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SUMMARY 

Following  increased  activity  in  underwater  ccoustics  research  during  World  War  II,  wave  theory 
modeling  of  acoustic  propagation  has  become  an  important  part  of  sonar  system  design  and  evalua¬ 
tion.  Since  then,  wave  propagation  analysis  of  sound  in  the  ocean  has  been  applied  to  the  interpreta¬ 
tion  of  marine  geophysical  measurements  and  is  now  an  established  tool  in  underwater  oil  explor¬ 
ation.  This  paper  prevents  a  survey  of  the  development  of  wave  theory  analysis  for  acoustic 
propagation  in  a  homogeneous  layer  of  fluid  overlying  a  homogeneous  fluid  half-space.  Although 
this  model  of  the  ocean  is  highly  simplified,  the  analysis  describes  the  important  features  of  the 
sound  field  and  forms  the  basis  for  constructing  more  sophisticated  models  of  the  ocean  environ¬ 
ment. 

Seqinning  with  the  axially  symmetric  wave  equation,  the  elementary  cylindrical  wave  solution 
is  found  by  separation  of  variables.  Following  the  technique  developed  by  Lamb  1,  through  an 
integral  summation  of  these  elementary  wave  functions  the  point  harmonic  source  is  replaced  by  a 
horizontal  boundary  plane.  In  this  manner,  the  solution  of  a  homogeneous  differential  equation 
rather  than  an  inhomogeneous  equation  is  required.  The  resulting  solution,  a  Hankel  transform 
integral,  may  be  evaluated  by  residue  theory  in  the  complex  wave  numl  er  plane. 

In  general,  the  form  cf  the  solution  depends  upon  the  path  of  the  branch  cut.  Pekeris  2 
employed  a  cut  that  yields  a  three  part  solution  in  which  each  part  corresponds  to  a  distinct  type  of 
propagation.  Thus,  a  finite  number  of  real  poles  gives  rise  to  a  set  of  modes  representing  waves 
trapped  in  the  layer.  An  infinite  number  of  modes  with  complex  eigenvalues  describe  energy  leakage 
out  of  the  layer  and,  finally,  a  branch  line  integral  contains  the  contribution  of  the  interface  wave. 
After  the  Pekeris  cut  has  been  examined,  the  solution  resulting  from  a  different  branch  cut  path 
proposed  by  Ewing,  Jardetsky  and  Press  3  is  developed.  This  cut,  the  EJP  branch  cut,  eliminates  the 
complex  modes;  however,  the  simplicity  of  the  previous  branch  line  integral  is  lost. 

Detailed  investigation  of  the  wave  that  travels  along  the  interface  between  the  layer  an d  the 
half-space  requires  a  slight  modification  of  the  model.  To  isolate  the  effects  of  this  wave,  the  surface 
of  the  fluid  layer  is  raised  to  infinity.  The  resulting  integral  is  asymptotically  evaluated  by  the 
method  of  steepest  descent  and  the  character  of  the  interface  wave  is  discussed.  In  addition,  these 
same  approximation  techniques  are  applied  to  the  Pekeris  branch  line  integral. 

Once  the  solution  techniques  for  the  layered  half-space  have  been  developed,  they  are  applied 
to  a  practical  shallow  ocean  propagation  problem.  Through  this  application,  the  significance  of  each 
type  of  propagation  is  illustrated. 

Following  the  analysis,  a  summary  of  the  types  cf  wave  propagation  in  this  simpie  ocean 
model  is  presented.  The  most  significant  propagation  component  at  long  range  consists  of  trapped 
waves  which  lose  no  energy  to  the  half-space  below.  At  short  range,  two  additional  .ypes  of 
propagation  are  evident.  One  of  these  field  components  is  the  result  of  travelling  waves  that  leak  out 
of  the  layer  thereby  decaying  rapidly  with  increasing  range.  The  other  short  range  component,  the 
interface  wave,  is  strongly  excited  only  when  a  mode  is  near  cutoff  in  the  layer. 
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DESCRIPTION  OF  THE  PHYSICAL  PROBLEM 

In  1948,  C.L.  Pekeris  published2  the  solution  to  the  probltm  of  the  oropegation  of  erouttic 
waves  in  a  fluid  layer  in  eontKt  with  a  fluid  half-space.  This  study  arose  as  a  consequence  of  a  series 
of  measure  men  a  performed  by  Worzel  end  Ewing*  employing  explosives  in  shallow  ocean  water. 
Ewing  noted  that  the  spectrum  of  the  energy  arriving  at  a  fixed  point  distent  from  the  source  was 
dispersed  or  spread  in  time.  In  other  words,  each  frequency  component  travelled  with  its  own 
characteristic  speed.  Using  the  theoretical  work  of  Gutenberg®  and  Stoneley®  on  the  dispersion  of 
Love  waves,  Pekeris  formulated  a  model  of  the  ocean  environment  corresponding  to  those  authors' 
work.  Since  Love?  had  demonstrated  the  existence  of  trapped,  dispersive  sheer  waves  in  an  elastic 
layer  over  an  elastic  half  space,  Pekeris  chose  to  analyze  the  equivalent  problem  for  two  fluids. 

The  resultant  analysis  successfully  explained  the  observed  dispersion  and,  subsequently,  this 
two  fluid  model  became  the  foundation  for  most  shallow  water  acrustic  modelling.  Inter  Improve¬ 
ments  in  the  analysis  included  allowing  the  sound  speed  to  vary  in  the  layer,  adding  more  layers  and 
introducing  elastic  properties  in  the  lower  half  space.  In  spite  of  these  extensions  to  the  theory,  the 
basic  solution  technique  remains  unchanged. 

In  this  paper,  we  shall  consider  the  analysis  of  the  two  fluid  model  in  detail.  Specifically,  the 
model  consists  of  one  homogeneous  fluid  layer  overlying  a  half-space  containing  a  second  homoge¬ 
neous  fluid.  As  a  consequence  of  the  fluid  media,  shear  waves  do  not  exist  in  this  model.  However, 
if  elastic  substances  were  substituted  for  the  fluids,  the  solution  for  the  shear  wave  field  would  be 
conceptually  identical  to  the  compressions!  wave  solution.  Excitation  of  the  field  is  by  means  of  a 
point  source  of  harmonic  waves  located  in  the  fluid  layer.  Also,  the  combination  of  a  point  source 
and  parallel,  planar  boundaries  leads  naturally  to  a  solution  expressed  in  cylindrical  coordinates. 

Normally  the  ocean  floor  sediment  possesses  a  higher  sound  speed  than  the  ocean  water  itself. 
Therefore,  in  our  analysis,  we  will  assume  that  the  sound  speed  of  the  fluid  in  the  half-space  is 
higher  than  the  sound  speed  in  the  layer.  While  there  are  instances  of  a  lower  sound  speed  in  the 
sediments  compared  to  the  sound  speed  in  the  overlying  water,  usually  only  the  upper,  highly 
saturated  sediment  layers  exhibit  this  behavior.  For  this  reason  we  will  confine  our  analysis  to  the 
problem  of  a  higher  sound  speed  in  the  half-space.  Additionally,  the  higher  sound  speed  in  the 
underlying  half-space  permits  the  existence  of  trapped  waves  which  dc.'Mnate  the  field  at  long 
range. 
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ANALYSIS 


Axially  Symmetric  Wave  Equation 

In  order  to  describe  a  point  source  in  a  medium  with  parallel,  planar  boundaries,  the  wave 
equation  can  be  expressed  in  cylindrical  coordinates.  Moreover,  since  the  analysis  does  not  contain 
any  azimuthal  variation,  the  axially  symmetric  form  of  the  wave  equation  is  sufficient,  i.e., 


where 

c  ■  sound  speed 
r  *  radial  distance  from  z  axis 
t  *  time 
z  -  depth 

i h  *  velocity  potential. 

Use  of  velocity  potential  reduces  the  equation  of  motion  to  a  scalar  rather  than  a  vector  equation. 


The  particle  velocity  is  related  to  the  velocity  potential  by. 


-V* 


3u  . 

TT* 


3v 

3t 


+  k 


dw 

w 


where 

D  «  displacement  vector 
u  *  x  component  of  displacement 
v  •  y  component  of  displacement 
w  *  z  component  of  displacement. 


(2) 


Furthermore,  the  volume  expansion  or  dilatation  of  a  small  element  of  fluid  is  given  by  the  di¬ 
vergence  of  the  displacement  and  is  related  to  the  pressure  in  the  fluid  by, 

P  -  -  E  V  •  D  (3) 


where 

P  =  pressure 

E  =  bulk  modulus  of  fluid. 


Also,  the  sound  speed  in  the  fluid  is  given  by, 


(4) 


where 

p  =  density. 
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Taking  the  divergence  of  equation  (2),  we  obtain, 

V2*  -  -  A  (v  .  o) 

Using  equation  (3)  for  the  pressure  and  equation  (4)  for  the  sound  speed  in  this  last  result,  we 
find, 


However,  equation  (1)  and  this  result  imply, 

P  ■  (6) 


Since  the  point  source  has  been  assumed  to  be  harmonic,  the  time  dependence  of  tf.e  potential 
function  can  be  written  as  follows, 

^  .  Aio)t  (r,  z)  (6) 


where 

i  ■  -s/  '-l 

cj,  angular  frequency  *  2vf 
<l>  *  harmonic  potential. 

Introducing  this  last  assumption  into  the  wave  equation,  i.e.,  equation  (1),  we  have. 


1  a  /  3<1>  \ 
r  dr  \  r~5r"  / 


92* 


CJ* 


+  Si?  +ir* 


(71 


The  assumed  time  variation  of  equation  (6)  also  leads  to  a  simpler  relationship  between  the  acoustic 
pressure  and  the  potential  function  than  indicated  by  equation  (5): 

p  »  iwp  4*  (8) 

Applying  separation  of  variables  to  the  time  independent  differential  equation  (7)  we  may  write 
a  range  equation, 

2  d2R  dR  ,  o  2d  n 
'  ^  +  'W+-,*'2*  -0 


and  a  depth  equation, 


where 

R(r)  =  function  of  the  radius,  r 
Z(z)  =  function  of  depth,  z 
7  =  horizontal  wave  number 


8 


NADC-81284-30 

The  range  equation  is  a  form  of  Bessel's  equation  and  the  solution,  bounded  as  the  range  increases 
to  infinity,  is  given  by, 

R=J0(7r)  (10) 

where  J 0{p)  ~  zero  order  Bessel  function  of  the  first  kind.  On  the  other  hand,  both  solutions  of 
the  depth  equation  (9)  must  be  retained  until  we  consider  the  fluid  layer  problem.  These  solutions 
are  of  the  form, 


2  =  e±  '0Z 


where 


p,  vertical  wave  number  = 


_  y2 
c2  7 


Therefore,  an  elementary  solution  the  axial  symmetric  wave  equation  can  be  written  by  sub¬ 
stituting  the  product  of  equations  (10)  and  (11)  for  the  harmonic  potential,  4>,  in  equation  (6). 
The  result  of  this  operation  is, 

<//  =  e'«t  JQ  (7r)  e^fr  (13) 

Observe  that  the  vertical  wave  number,  p,  is  defined  in  equation  (12)  by  a  square  root  and,  there¬ 
fore,  can  be  either  positive  or  negative.  In  order  to  select  the  proper  sign,  we  consider  only  the 
exponential  factors  of  equation  (13)  which  can  be  combined  as  follows, 


* 

J0<7r) 


=  ei(cut 1  & 


If  the  wave  number,  (3,  is  real,  the  negative  sign  indicates  a  wave  travelling  in  the  positive  z  direc¬ 
tion  while  the  positive  sign  refers  to  a  wave  travelling  in  the  negative  z  direction.  In  some  regions 
of  the  medium,  waves  travelling  in  both  directions  are  physically  reasonable;  however,  in  a  region 
beyond  any  reflectors,  there  is  nothing  to  cause  the  wave  to  return  toward  the  source.  For  this 
case,  the  sign  of  the  wave  number  must  be  selected  to  guarantee  waves  which  travel  away  from 
the  source.  In  the  event  that  oj2/c2  is  less  than  y2,  equation  (12)  produces  an  imaginary  value  for 
j3  and  the  depth  function,  Z,  is  exponential  in  behavior.  Thus,  the  sign  of  the  exponent  of  equa¬ 
tion  (14)  must  he  chosen  to  insure  a  physically  reasonable  decrease  in  the  amplitude  as  we  depart 
from  the  source.  This  consideration  will  be  discussed  in  greater  detail  in  connection  with  the  im¬ 
proper  or  leaky  modes  of  vibration. 

The  technique  we  shail  employ  to  develop  the  solution  to  the  layered  problem  is  generally 
attributed  to  Lamb^,  A  point  source  generates  a  spherical  wave,  at  least  before  any  boundary 
interaction  occurs,  while  a  cylindrical  coordinate  system  is  convenient  for  the  description  of  the 
boundary  conditions.  Therefore,  we  shall  expand  a  spherical  wave  in  terms  of  the  elementary  cylin¬ 
drical  wave  functions  given  by  equation  (13).  This  may  be  accomplished  by  the  application  of  the 
Hankel  transform  accord:  ,g  to  the  procedure  outlined  in  Appendix  A.  The  result  of  this  develop¬ 
ment,  equation  (A-12),  i- 


=  piout 


J0  (7r)  e 

iF 


rr  tr  arr';?-  * 
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Initially,  we  shall  consider  the  source  to  be  located  in  the  plane,  2  =  0.  Following  Pekeris'  tech¬ 
niques,  we  represent  the  source  as  a  set  of  boundary  conditions  rather  than  as  an  inhomogeneous 
source  term  in  the  original  differential  equation.  Thus,  we  require  continuity  of  the  pressure  and  a 
discontinuity  in  the  particle  velocity  across  the  plane,  z  =  0,  such  that  the  depth  function  assumes 
the  form, 

2  =  Ze±i/3z 
*  i/3 6 

This  representation  of  the  depth  function  guarantees  that  the  z  dependence  of  the  resulting  po¬ 
tential  function  is  identical  to  the  integral  indicated  in  equation  (15).  Therefore,  the  general  form 
for  the  depth  equation  is, 

2+  =  Ct  e-'frfzX)) 

Z_  =  C2  e$z  (z  <  0) 

Note  that  the  signs  in  the  exponents  have  been  selected  to  insure  that  the  resulting  waves  diverge 
from  the  source.  Continuity  of  pressure  along  the  plane,  z  =  0,  yields, 

C1  =  c2  =  jjy  (16) 

Further,  the  discontinuity  in  the  vertical  component  of  particle  velocity  can  be  expressed,  using 
equation  (2),  in  the  form, 

dZ_  dZ+ 

-ar  -  d t  ,,7) 

o  o 

Therefore,  we  shall  introduce  a  point  source  into  our  analysis  by  placing  a  horizontal  boundary 
plane  through  the  source  and  requiring  that  the  two  conditions  expressed  by  equations  (16)  and 
(17)  be  satisfied  on  this  source  plane. 


Formulation  of  the  Boundary  Value  Problem 


Since  we  have  replaced  the  point  source  by  a  boundary  surface,  we  can  formulate  the  prob¬ 
lem  of  a  fluid  layer  overlying  a  fluid  half- space*.  We  shall  assume  that  the  fluid  layer  is  bounded 
above  (see  figure  1)  by  air  and  below  by  a  second  fluid  with  a  higher  sound  speed  than  the  layer. 
Since  the  density  of  air  is  so  much  lower  than  th3t  of  water,  the  fluid  in  the  layer,  we  shall  assume 
that  the  pressure  on  the  upper  surface  of  the  layer,  i.e.,  the  plane,  z  =  0,  is  equal  to  zero.  Because 
the  sound  speed  in  the  layer  is  not  equal  to  the  sound  speed  in  the  underlying  half-space,  the 
vertical  wave  number,  (3,  will  have  a  different  value  in  each  region.  The  expressions  for  the  wave 
numbers  can  be  written  from  equation  ( 1 2)  as, 


10 


FIGURE  1.  -  Geometry  of  the  layered  half-space 
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where 

c-j  =  sound  speed  in  the  layer 
C2  =  sound  speed  in  the  half-space. 

In  addition,  the  depth  functions,  equation  (11),  may  be  expressed  in  the  form. 

Region  I  :  Zy  =  C-j  sin  /?iz  +  C2cos0-|z 

Region  II  :  Z2  =  C3sin0iz  +•  C4COS/31Z  (19) 

Region  III:  Z3  =  Cg  e  '^2Z  +  Cg  e '^2 
where  the  regions  are  defined  in  figure  1. 

Recalling  the  previous  discussion  concerning  the  sign  of  the  exponent  in  equation  (14),  we  ob¬ 
serve  that  the  second  term  in  the  expression  for  Z3  represents  waves  converging  on  the  source  from 
the  half-space.  Nothing  in  the  half-space  can  produce  waves  travelling  toward  the  source,  therefore, 
this  second  term  is  meaningless,  or 

C6  =  0  (20) 

However,  if  7 2  is  greater  than  ,  then,  by  equation  (18),  fa  is  imaginary  and  must  be  writ¬ 

ten  as, 

fa  -  (21) 

Notice  that  the  negative  square  root  has  been  selected  in  order  that  when  this  last  result  is  sub¬ 
stituted  into  the  third  depth  function  expression,  equation  (19),  the  first  term  amplitude  decreases 
exponentially  with  increasing  depth  in  the  half-space. 

In  view  of  the  previous  remarks  concerning  the  air-water  interface  and  the  representation  of 
the  source  as  a  boundary,  we  may  express  the  conditions  along  the  horizontal  boundaries  of  the 
three  regions  as  follows: 

1)  at  z  =  0  the  pressure  is  equal  to  zero,  or, 

Z^OJ-O;  (22) 

2)  at  the  source  depth,  z  =  d,  in  accordance  with  equations  (16)  and  (17),  we  have, 

Z-|  (d)  -  Z2(d) 

for  continuity  of  pressure,  and, 

dZi  (d)  dZ9  (d) 

"dl-  -  ~a—  =  37 


fo"  the  discontinuity  in  particle  velocity;  and, 
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3)  at  the  interface,  Z  =  H,  separating  the  layer  from  the  half-spiace,  we  express  the  pressure 
by  means  of  equation  <8), 

Pi  Z2  (H)  =  P2  Z3  (H) 

for  continuity  of  pressure,  and, 

dZ2  (H)  dZ3(H) 
dz  ~  dz 

for  continuity  of  particle  velocity. 


Introducing  the  depth  functions  given  by  equations  (19)  into  the  above  boundary  conditions,  we 
obtain  four  equations  in  four  unknowns.  Since  the  first  condition,  equation  (22),  leads  immedi¬ 
ately  to  C2  =  0,  we  may  write, 

Ci  sin  0i  d  -  C3  sin  0i  d  -  C4  cos  0i  d  =0 

01  Ci  cos  0i  d  -  0i  C3  cos  0i  d  +  0i  C4  sin  0i  d  =  2y 

Pi  C3  sin  0i  H  +  pi  C4  cos  0i  H  -  P2  Cg  e  '^2^  =  q 

01  C3  cos  0i  H  -  0i  C4  sin  0i  h  +  i02  C5  e  '^2^  =  q 


The  constants  are  evaluated  by  solving  this  system  of  equations  simultaneously.  Thus, 

27  [~~Pl  cos  [01  (H-d)]  +  ib02  sin  [0i  (H-d)] 

^  01  01  cos  0i  H  +  ib0p  sin  0i  H 


c3  - 

C4  = 

C5  = 


27  sin  0i  d 

37 

27  sin  0i  d 


01  iiii  0i  H  -  ib02  cos  0i  H 
01  cos  0i  H  +  ib02  sin  0i  H 


27be'^  ^  sin  0i  d 
01  cos  0i  H  +  ib02  sin  0i  H 


where 


b  =  pi/p2 


By  substituting  the  values  of  the  above  constants  in  equations  (19)  and  constructing  potential 
functions  in  the  form  of  equation  (13),  the  velocity  potential  for  each  region  can  be  written  as, 


^1  ~  2eiojt  f  JQ  (7r)  7 

J a 


sin  0iz 
01 


PI  cos  pin  +  1DP2  sin  pi 


d7 


(0  <  z  <  d) 


(23) 
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^2 


»  f  sin0id 

2e»«t  J  Jcl  yr)y 

n 


sin  0id  |  (3i  cos  t/3-| ( H-z) ]  +  ib02  sin  (H-z)] 


cos  0i  H  +  ib02  sin  0i>H 


d7 


(d  <  2  <  H) 


(24) 


</'3 


f 


2e'wt  j  j0(7r)7 
o 

(H  <z) 


b  sin  0ide 


;02<2-H) 


— I 


01  cos  0i  H  +  ib02  sin  0i  H 


d7 


(25) 


To  complete  the  analysis,  we  must  evaluate  the  above  integrals  and,  for  this  purpose,  we  will  trans¬ 
form  them  into  complex  contour  integrals.  We  observe,  however,  that  the  expression  for  \jj 2  is  iden¬ 
tical  to  the  expression  for  4/-\  if  z  and  d  are  interchanged.  Thus,  in  order  to  compute  the  field 
anywhere  in  the  fluid  layer,  we  need  only  evaluate  equation  (23)  or  equation  (24).  These  integrals 
can  be  evaluated  numerically  when  r  is  small  and  the  wavelength  is  comparable  to  the  layer  depth, 
H.  In  all  other  cases,  the  integrand  oscillates  rapidly  and  is  difficult  to  integrate  numerically. 

Integration  in  the  Complex  Wave  Number  Plane 


Before  integrating  equations  (23),  (24),  and  (25),  it  is  useful  to  replace  the  Bessel  function, 
J0,  bv  Hankel  functions^.  The  Hankel  functions  are  related  to  the  Bessel  function  as  follows. 


J0<7r)  -  y  ^H0<1>(7r)  +  H0<2>(7r)J 

(26) 

Hq^M-  -H0<2){-7r) 

(27) 

where 

H0^  -  Hankel  function  of  the  first  kind 

H0^2)  «  Hankel  function  of  the  second  kind 

For  large  values  of  the  argument,  7r,  the  Hankel  functions  may  be  replaced  by  asymptotic  approxi- 

mations  of  the  form, 

H0n,trl  *  ^e'  ^r"7) 

(28) 

h0I2,m  = 

(29) 

where 


7r  »  1 
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First,  let  us  rewrite  equation  (24)  through  the  use  of  equation  (26), 


=  e 


lOJt 


/  H0^(7r)y  G(/3-j,  02>  d7 
o 


where 


OO 

f  H0'2^(7r)7  G(]3i,  |32>  d7 


>  (d  <  z  <  H) 


sin  0i  d  pi  cos  (0i(H-z)]  +  ib02  sin  [/3-| (H-z) ]  | 
G(01- 02^  I _  01  cos  0i  H  +  ib02  sin  0i  H 

Relation  (27)  can  be  employed  to  rewrite  this  integral  as  follows, 


(30) 


OO 

'P 2  =  e'wt  f  H0^2^(7r)7  G(01,  02>  d7 


(31] 


In  order  to  establish  the  behavior  of  the  integrand  in  the  7-plane,  we  recall  that  the  defini¬ 
tions  of  0i  and  02  from  equations  (18)  lead  to  the  branch  points  given  below, 


(32) 


Further,  the  poles  of  the  integrand  are  located  at  the  zeros  of  the  denominator  of  equation  (30), 
or, 


01  cos  0i  H  +  ib02  sin  0i  H  =  0 


(33) 


Wo  will  show  that  the  pole  at  0i  =0  lies  outside  the  contour  used  to  evaluate  equation  (31),  thus, 
it  contributes  nothing  to  the  solution. 


Pekeris^  proposed  the  contour  illustrated  in  figure  2,  in  which  the  branch  cuts  extend  down¬ 
ward  from  the  branch  points.  Observe  that  as  the  radius  of  the  circular  arc  contours,  Ci  through 
C4,  becomes  infinite,  the  path,  A,  becomes  the  integration  path  for  the  desired  potential  integral, 
equation  (31).  We  express  this  relationship  as  follows, 


i>2 


^2 

eiuJt 


A 


/•  /• 

c1234  B] 


( residues)  n 
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Pekeris  Cut 
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In  the  above  expression  the  sum  includes  the  residues  of  all  the  poles  enclosed  by  the  contour. 
Since  the  Henkel  function  decays  exponentially  for  large  magnitudes  of  7  in  the  third  and  fourth 
quadrants  as  described  by  equation  (29),  the  integrals  over  the  paths  C|,  C2.  C3  and  C4  approach 
zero  for  large  radius.  (A  formal  proof  for  this  behavior  is  presented  in  Appendix  B.)  Also,  from 
equation  (30),  we  notice  that  G(0i,  02 )  is  an  even  function  of  0i,  hence,  the  integral  along  one  side 
of  the  branch  cut  cancels  that  along  the  other  side.  Thus,  the  solution  reduces  to, 


<l>2  a  -2ir\  2  residues  +  /  (34) 

h 

The  residues  of  equation  (31)  can  be  found  by  rewriting  ths  intcgrond  as, 


where 


D(y)  =  denominator  of  equation  (30)  without  the  (3y  factor 
N(-y)  =  remainder  of  the  integrand. 


Therefore, 


2iri(  residue)  n  =  2rri  ^  ^0(7^^ 


2rri 


=  2iri 


7=7n 

N(7)  +  (7-7n>3^  N(y) 
N(7n» 


7=7n 


97 


0(7) 


7=7n 


By  performing  the  indicated  operations,  we  find, 

(2) _ fil  sin  /^dsinfliz 


( residue)  n  =  H0uM7nr) 


H0i  -  cos  0i  H  sin  0i H  -  sin^  0i  H  tan  0i H 


(35) 


where 

01  is  evaluated  for  7  =  7n 

7n  is  the  value  of  7  at  the  nth  pole. 
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Furthermore,  we  observe  thet  02  has  been  eliminated  from  this  equation  by  using  equation  (33). 

By  replacing  the  Hankei  function  in  equation  (35)  by  its  asvmptotic  form,  equation  (29), 
the  physical  significance  of  each  residue  or  normal  mode  cm  be  determined.  The  range  and  time 
dependent  factors  of  each  mode  may  be  combined  so  thet  the  mode  it  asymptotically  proportional 
to  the  following  expression, 


( residue) n  <*  — - j~  e',CJt  ~Tnr) 


(36) 


where 


CJ 


(37) 


Hence,  each  mode  represents  a  travelling  cylindrical  wave  with  a  phase  velocity,  cn.  In  the  event 
that  7n  is  real,  the  wave  suffers  no  attenuation,  only  the  cylindrical  spreading  given  by  the  r-/* 
factor.  In  other  words,  the  wave  is  completely  trappea  in  the  fluid  layer.  On  the  other  hand,  if  yn 
is  complex,  the  wave  will  be  damped  because  of  the  appearance  of  a  real  exponential  factor  in 
equation  (36).  We  obse'-ve  that  the  contour  in  figure  2  was  chosen  so  that  the  imaginary  part  cf 
7  will  always  be  negative,  thus  insuring  that  a  real  exponential  factor  in  equation  (36),  if  it  exists, 
decays  with  increasing  range.  This  explains  the  selection  of  H0^)  in  equation  (31)  since  this  func¬ 
tion  passes  to  zero  with  large  negative,  imaginary  argument. 

Proper  evaluation  of  the  uranch  line  integral,  depends  on  the  choice  of  the  sign  of  02  on  «*»ch 
side  of  th  „•  branch  cut.  txpressing  the  horizontal  wave  number,  y,  as  a  complex  number, 

7  »  7r  +  '  7j 


where 

7r  =  real  part  of  y 
7j  =  imaginary  part  of  y, 
we  may  represent  02  employing  equation  (18), 

2 

022  -  "Tr 2  +  7i2  -  2i  vr7i  (38) 

c2^ 

Along  the  branch  cut,  7r  =  u;'C;2,  and  7j  ranges  from  0  to  therefore, 

P22  :  ‘“i2  ~  2’7j  ~  (-  «  <  7j  <  0) 

'*2 

This  parabola  is  shown  as  contour  B,  in  figure  3  jnd  A+  and  A_  represent  the  two  roots  of  022< 
i.e.. 


A-  -  - 
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As  demonstrated  in  Appendix  C,  the  real  eigenvalues,  7n,  corresponding  to  trapped  inodes  are  all 
located  on  the  real  axis  of  the  7-plane  to  the  right  of  w/c2.  Substituting  7j  *  0  and  yr  >  w/c2 
into  equation  (38),  we  find  that  this  portion  of  the  real  axis  in  the  7-plane  transforms  to  either 
the  positive  or  negative  imaginary  axis  of  the  02*plane  since, 

(j22  .  _  yr2  <  0  (39) 

c2* 

If  we  recall  the  depth  function  in  the  half-space  below  the  layer,  i.e.,  equations  (19)  and  (20), 

23  *  C5e'^22  (40) 


we  conclude  that,  if  02  is  imaginary,  it  must  be  negative  in  order  to  guarantee  the  physically  reason¬ 
able  exponential  decay  with  increasing  depth  in  the  half-space.  Therefore,  we  select  the  region  be¬ 
low  contours  A.  and  A+  in  figure  3  as  the  branch  for  02 .  By  substituting  the  appropriate  values 
for  7r  and  7;  into  equation  (38),  we  can  show  that  the  path  In  the  7-plane,  labelled  number  2  in 
figure  3,  transforms  to  the  positive  real  axis  of  the  02  plane.  The  transformation  between  the 
7-plsne  and  the  02-plane  is  illustrated  in  figure  3. 

On  the  basis  of  figure  3,  we  can  summarize  the  behavior  of  02  in  the  7  plane  as  follows: 


Region  1:  02r  >  0,  02i  0 

(use  positive  root  for  02) 

(41) 

Region  1 1 :  02r  <  0*  02j  <  0 

(use  negative  root  for  /?2> 

(42) 

Region  III:  02r  >  0-  02j  <  0 

(43) 

where 

02r  »  real  part  of  02 
02i  *  imaginary  part  of  02 


Now  that  we  have  determined  which  value  of  02  to  select  on  each  side  of  the  branch  cut,  the 
branch  line  integral  can  be  specified.  As  shown  in  figure  2,  the  psth  passes  up  along  the  left  side 
of  the  cut  and  down  along  the  right,  so  that  the  integral  in  equation  (34)  becomes, 
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By  substituting  the  expression  for  G  (0j,  02),  equation  (30),  into  this  last  result  and  simplifying,  the 
branch  line  integral  can  be  written  as, 


02  sin  0id  sin  (Jy  z 


0j2cos20iH  +  b2  022  sin20jH 


dy 

I 


(44) 


Therefore,  the  complete  solution,  equation  (34),  is  the  combination  of  this  branch  line  integral  and 
the  sum  of  modes  given  individually  by  equation  (35).  Hence, 


4»(r,z)  ■  -2»iT3  H0^  (ynr)  - — - — - 1'  ■ 

"  H0i  -  cos 0i H  sin  0^H  -  b2  sin  0-|H  tan  0^H 


0y  sin  0jd  sin  0j  z 


<jj 

c2 


-2ib 


r  2  r 

J  H0(2)(7r)7  - 

i«o  L 

«2 


02  sin  0^d  sin  0^  z 


0y 2  cos2  JyH  +  b2  022  sin2  0yH 


r] 


d7  (45) 


According  to  Appendix  C,  the  roots  of  the  characteristic  equation,  equation  (33),  either  lie 
along  the  real  axis  between  w/c2  *nd  gj/c]  or  in  region  I  of  the  7-plane  shown  in  figure  3.  The 
former  set  of  solutions  is  a  finite  set  and  represents  waves  with  energy  trapped  in  the  layer  while 
the  latter  set  of  solutions  is  infinite.  Each  mode  of  this  infinite  set  spreads  cylindrically,  i.e.,  its 
amplitude  varies  as  r1/4,  as  do  the  trapped  modes.  However,  these  modes,  which  are  not  trapped, 
also  decay  exponentially  with  increasing  range  because  the  eigenvalue,  7n,  has  a  negative  imaginary 
part.  These  modes  corresponding  to  complex  eigenvalues  are  called  “leaky"  or  "improper"  modes. 
Such  modes  are  leaky  because  the  energy  is  not  completely  trapped  in  the  layer  and  improper 
because  of  their  behavior  as  depth  in  the  half-space  increases. 

Consider  the  vertical  wave  number,  02-  'n  region  I  of  the  7-plane,  the  region  in  which  these 
improper  modes  are  found.  In  this  region,  both  the  real  and  imaginary  parts  of  02  are  positive  as 
expressed  by  equation  (41).  Hence,  the  depth  function  in  the  half-space,  i.e.,  equation  (40),  in 
combination  with  the  time  dependence  is  of  the  form, 

Z3  =  eiwt  e'^2*  *  ei(^2r*> 

This  last  expression  describes  a  wave  diverging  from  the  source  as  it  should;  however,  the  last 
factor  indicates  that  it  increases  exponentially  with  depth.  Such  unphysical  behavior  for  a  mode 
leads  to  its  being  termed  an  improper  mode.  It  has  been  found  that  the  sum  of  these  improper 
modes  is  well  behaved  at  large  depths.  Indeed,  useful  results  have  been  obtained  for  the  field  in 
the  layer  by  adding  these  modes  and  the  trapped  modes^-  11. 

As  we  have  observed,  the  improper  modes  are  found  to  lie  in  region  I  of  the  7-pl8ne  and  ex¬ 
amination  of  figure  3  suggests  that  this  region  could  be  eliminated  by  re-routing  the  branch  cut. 
If  the  cut  is  taken  from  w/c2  along  the  real  axis  to  the  origin  and  then  down  the  negative  imaginary 
axis  to  infinity,  the  contour,  A+,  in  the  02*plane  would  lie  along  the  positive  real  axis  and  A. 
would  coincide  with  the  negative  real  axis.  As  a  result,  region  I  disappears  and  the  only  modes 
remaining  corresoond  to  trapped  modes  with  poles  on  the  real  axis. 
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This  cut  was  proposed  by  Ewing,  Jardettky  and  Press*  and  is  called  the  EJP  cut  in  this  paper 
while  the  cut  considered  previously  is  called  the  Pekeris  cut.  The  branch  line  integral  is  identical 
in  form  to  the  Pekeris  branch  line  integral,  equation  (44);  however,  the  limits  of  integration  are 
modified  to  follow  the  EJP  cut.  Thus,  the  EJP  branch  line  integral  is  given  by, 


-2ib 


Jd. 

c2 

f  H0^  (*yr)  y 


*1 


<32  sin  ftid  sin  (S]  z 


*  cos*  (iy  H  ♦  b*  02*  sin*  fiy 


H 


d7 


(46) 


The  EJP  cut  is  shown  as  path  2  in  figure  3.  Although  the  prospect  of  not  having  to  locate  and 
evaluate  complex  eigenvalues  is  attractive,  the  branch  line  integral  for  the  EJP  cut  is  much  more 
difficult  to  evaluate  numerically  than  the  Pekeris  branch  tine  integral,  particularly  at  high  fre¬ 
quencies1*.  This  is  true  because  the  Henkel  function  oscillates  much  more  rapidly  when  the  real 
part  of  its  argument  is  varied,  as  in  the  final  segment  of  the  EJP  cut  along  the  real  axis,  than  when 
the  imaginary  part  is  changed. 


Reflection  of  a  Spherical  Wave 


We  have  completed  the  solution  to  the  problem  of  a  layer  overlying  a  half-space,  but  we 
have  not  interpreted  the  branch  line  integral,  equation  (44),  physically.  In  order  to  gain  some 
insight  into  the  interaction  of  waves  with  the  interface,  it  is  necessary  to  simplify  the  problem 
somewhat  To  accomplish  this  purpose,  we  shall  let  the  layer  thickness  become  infinite.  In  this 
case,  the  problem  reduces  to  one  of  two  half-spaces  and.  as  a  result,  the  complications  of  waves 
returning  from  surface  reflections  are  avoided. 


Anothv..  benefit  of  this  simplification  is  that  the  lower  half-space  need  no  longer  be  homogene¬ 
ous  provided  that  it  can  be  characterized  by  a  plane  wave  reflection  coefficient.  We  gain  this  gen¬ 
eralization  by  expanding  the  spherical  wave  from  a  point  source  in  terms  of  plane  waves  rather 
than  in  terms  of  cylindrical  waves.  From  the  derivation  in  Appendix  D.  a  spherical  wave  can  be 
expressed  as  equation  (D-6)  which  is  repeated  below. 


e~ikR  _1_ 
R  *  2ir 


(47) 


where 

K,  vector  wave  number 
R,  slant  range 
<ji,  vector  range 
x,  y.  z 
7x*  7y.  J 


*  7x'  +  7yj  ♦  <3k 

-  /x*  +  y*  +  z* 

=  xi  +  yj  +  zk 

-  position  coordinates  of  a  point 
=  plane  wave  number  components 
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In  equation  (47),  the  factor  e-'^  '  ^  denotes  a  plane  wave  with  the  direction  of  the  normal  to 
the  wave  front  given  by  the  K  vector.  Also,  since  the  source  is  at  the  origin,  the  plane  wave  has 
travelled  through  the  distance  given  by  fl.  For  a  spherical  wave,  the  radial  direction  vector  is  always 
perpendicular  to  the  wave  front,  i.e.,  parallel  to  K,  hence,  the  dot  product  is  unnecessary  in  the 
exponent  of  the  left  hand  side  of  equation  (47). 

Employing  the  geometry  of  figure  4,  the  field  at  P  can  be  written  in  terms  of  the  direct  and 
reflected  paths, 


4V 


oo 

.-L 

2*  JJ  i if 


«1 

j/3  d7x  dTy 


f  fe'iK  •  (<R2a  + 

2 it  J  J  i/3 


(48) 


«2b) 


V(a)  d7x  d7y 


where 


V(a)  =  plane  wave  reflection  coefficient 
a  =  gra2ing  angle 

The  reflection  coefficient  has  been  introduced  in  the  second  term  of  equation  (48)  in  order  to 
account  for  the  reflection  of  the  plane  wave  at  the  boundary  between  the  half-spaces.  By  expand¬ 
ing  the  dot  products,  we  can  rewrite  the  two  terms  of  equation  (48)  as  follows. 


Y)  e-ij3(z-d) 
i/3 


d7x  d7v 


<j>2  = 


oo 

■iff 


j-i(7xX  +  7yY)  e_i,3(z+d) 
i/3 


V(a)  d7x  d7y 


(49) 


(50) 


where 

=  potential  for  the  direct  path 
<t>2  =  potential  for  the  reflected  path 

Ingard3  3  used  this  representation  to  summarize  Weyl's  work  on  the  reflectivity  of  spherical  waves; 
however,  he  adopted  the  angle  of  incidence  to  express  the  reflection  coefficient  and,  subsequently, 
to  transform  the  integral.  In  the  work  to  follow,  we  will  adopt  the  more  common  convention  in 
underwater  acoustics  and  employ  the  grazing  angle,  a. 

In  spite  of  the  appearance  of  figure  4,  the  function  <I>2  includes  the  entire  field  resulting  from 
the  interaction  of  t.ie  incident  wave  with  the  lower  half-space  —  not  only  the  simple  reflected  wave. 
Inhomogeneous  waves,  for  which  (3-]  is  imaginary,  are  discussed  in  Appendix  D  and  it  is  shown 
there  that  they  are  boundary  waves  travelling  along  the  interface  between  the  media.  These  waves 
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FIGURE  4.  -  Geometry  of  the  two  half-space  problem.  Paths  shown  are  normals  to  the 
wavefronts  of  a  direct  ( R i )  and  a  reflected  (R2)  plane  wave 
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are  shown  below  to  be  included  in  equation  (50).  In  order  to  investigate  the  behavior  of  this  com¬ 
posite  field,  we  will  transform  equation  (50)  by  employing  the  geometry  of  the  wave  number 
vector.  Figure  5  illustrates  the  relationship  between  the  wave  number,  k,  and  the  horizontal  and 
vertical  wave  numbers,  y  and  0,  and  the  grazing  angle  a.  Observe  in  equation  (47)  that  k  describes 
the  basic  spherical  wave  we  are  attempting  to  reproduce  and,  therefore,  is  a  constant.  Thus,  the 
extension  of  y  to  infinity  implies  that  j 3  and  a  must  become  complex. 

Following  the  same  procedure  illustrated  in  Appendix  D  to  obtain  equation  (D-7),  we  trans¬ 
form  equation  (50)  to  an  integration  over  y, 


$2 


e-i0l(z+d) 


V(a) y  d7 


(51) 


In  Appendix  E,  the  solution  for  a  homogeneous  lower  half-space  is  derived  by  solving  the 
corresponding  boundary  value  problem.  If  we  compare  equation  (E-4)  to  equation  (51),  we  ob¬ 
serve  that  the  two  results  are  equivalent  if  the  reflection  coefficient  is  defined  as, 


01  - 

01  +  b02 


(52) 


This  last  result  is,  in  fact,  the  plane  wave  reflection  coefficient  for  the  problem  of  two  homogene¬ 
ous  half-spaces.  A  more  familiar  form  of  this  result  can  be  written  in  terms  of  the  grazing  angle  a 
if  we  observe  that,  from  figure  5, 


7  =  k  cos  a  -  —  cos  a 
01  =  k  sin  a  -  sin  ot 


(53) 


Introducing  these  last  results  and  the  definitions  of  0i  and  02  Given  in  equations  (18)  into  equa¬ 
tion  (52),  the  reflection  coefficient  expression  becomes, 


V(a)  = 


sin 


sin 


Vci2 

— =:  -  cos^  a 
_ £21 _ 

a  +  b  aJ 


(54) 


£T 

c2: 


-  cos^  a 


This  last  result  is  equivalent  to  that  given  by,  for  example,  Clay  and  Medwin^. 

Equation  (51)  can  now  be  transformed  to  the  cr-plane  (the  complex  grazing  angle  plane) 
by  writing  the  equivalent  differential  element  from  equation  (53), 


d7  =  -  k  sin  a  d  a 


(55) 
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To  locate  the  path  of  integration  in  the  a-plane,  we  write, 

7  =  k  cos  (a1  +  i  a") 

(56) 

*  k  (cos  a' cosh  a”  -  i  sin  a  sinh  a"] 

where 

a  -  real  part  of  a 
a"  «  imaginary  part  of  a 
Thus,  it  may  be  shown  that, 

7  =  -“for  a'  =  jt  and  a"  -  ±  00 
7  =  0  for  a'  *  tt/2  and  a"  =  0 
7  =  +<»fora'  =  0  and  o"  =  ±  °° 

In  oraer  to  insure  that  there  is  an  exponential  decay  in  the  wave  with  increasing  depth  below  the 
source,  the  imaginary  part  of  0i  must  be  negative;  thus,  from  equations  (53)  we  write, 

01  =  k  sin  (a  +  ia") 

(57> 

=  k  [sin  a  cosh  a"  +  i  cos  a'  sinh  a." ] 

Therefore,  when  a  -  it.  a"  must  be  positive  and,  if  a'  *  0,  a"  must  be  negative.  Consequently, 
the  transformed  path  of  integration  starts  at  a  =  it  +  i  “,  moves  down  to  a  =  7r,  along  the  real 
axis  to  a  -  0,  and  ends  at  a  -■  -  i  This  path  is  the  reverse  of  path  C  in  figure  6. 

Applying  these  results  to  equation  (51),  the  field  integral  becomes, 

-j  oo 

<t>2  =  -y-  J'  H0^  (kr  cos  a)  e~ik  (z+c*)  sin  01  V(a)  cos  a  d  a 

7T  +  i  00 

where 

H0^)  =  Hankel  function  of  the  second  kind. 

Since  the  limits  of  integration  can  be  interchanged,  we  obtain, 


(kr  cos  ot)  e~ik(z+d)  s*n  a 


c 


V(a)  ccs  a  d  a 


where 

C  =  contour  shown  in  figure  6. 


(58) 
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FIGURE  6.  —  Configuration  of  the  complex  grazing  angle  plane  depicting  the  original 
integration  contour  C,  the  path  of  steepest  descent,  S,  through  the  saddle  point 
a0,  and  the  two  branch  lines,  B  and  B~ 
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Asymptotic  Expression  for  the  Reflected  Field 

Equation  (58)  is  the  exact  solution  for  the  reflected  field  of  figure  4;  however,  numerical 
evaluation  of  this  expression  is  difficult  because  the  integrand  may  oscillate  very  rapidly.  Hence, 
we  will  consider  an  asymptotic  approximation  for  the  solution  valid  for  large  values  of  kr.  This 
approximation  is  also  useful  for  providing  insight  into  the  physical  significance  of  the  reflected 
field.  In  the  following  analysis,  we  will  use  procedures  developed  for  electromagnetic  wave  piopa- 
gation^-  16(  although  our  solution  will  be  written  in  terms  of  the  grazing  angle  rather  than  the 
angle  of  incidence. 


First,  the  Hankel  function  is  replaced  by  its  asymptotic  approximation  and,  in  this  case,  we 
include  the  second  term  in  the  series  expansion^,  hence, 


Hq<2>  (kr  cos  a) 


-  J  —J— .  /,_ _ ] _ \ 

f  jrkrcosa  \  8ikrcosa  ) 


-i(kr  cos  a  -■?•) 


(59) 


Combining  the  exponential  factor  indicated  above  with  the  exponential  in  equation  (58),  we  ob¬ 
tain  an  exponential  of  the  form, 


«-ik  [r  cos  a  +  (z  +  d)  sin  a] 


If  we  designate  the  length  of  path  2  in  figure  4,  i.e.,  the  sum  of  R2a  and  Rjb.  by  R  and  let  the 
grazing  angle  shown  be  a0,  the  above  factor  becomes, 

e-ikR  cos  {ol-ol0)  (60) 

Introducing  these  changes  into  equation  (58),  we  obtain, 

*2  ’  e  ’ 4  /( 1  -  likrVoT^)  <r'kR  C°S(0"a°)  Vta)  do  (61 ) 

c 

Next,  we  examine  the  characteristics  of  the  integrand  in  the  complex  a-plane.  There  are  three 
singularities  of  importance: 


a)  the  two  branch  points  of  V(a),  or,  equation  (54),  given  by, 

^ J  — ^  -  cos2  a  =  0 

b)  the  poles  of  the  function  V(a),  and, 

c)  the  branch  point  located  at 

y/  cos  a  =  0 

Initially,  let  us  consider  the  branch  points  of  V(a), 


(62) 
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These  branch  points  and  their  associated  branch  cuts  can  be  treated  in  the  manner  employed  in 
the  layered  half-space  problem.  Referring  to  figure  3,  we  use  the  path  equivalent  of  the  EJP  cut, 
i.e„  the  positive  real  axis  in  the  Transformation  of  the  EJP  cut  from  the  7-plane  to  the 

a-plane  yields  path  B  in  figure  6.  Since  there  are  two  roots,  agp  and  ofQP-  t0  equation  (62),  a  sec¬ 
ond  branch  cut  must  exist  and  it  is  shown  as  B~  in  figure  6.  Further,  the  branch  is  selected  so  that 
the  imaginary  part  of  c-]^00  's  negative;  thus,  the  field  decays  with  depth  below  the  interface. 


The  poles  of  the  function  V(ot)  are  determined  by  the  zeroes  of  the  denominator  of  equa¬ 
tion  (54),  or, 


sin^  a  =  b^ 


(63' 


where 

a  ■  o'  +  ia" 

a'  =  real  part  of  a 

a"  -  imaginary  part  of  a 

There  are  three  roots  to  this  last  equation:  two  real  roots,  c*i  and  c*2>  and  one  imaginary  root, 
03,  as  illustrated  in  figure  6.  All  of  the  poles  resulting  from  these  roots  lie  outside  of  the  region 
defined  by  the  contour  integration;  therefore,  they  can  be  ignored.  Also,  the  branch  point  of 
v/~cos  a  located  at  a  =  ir/2  cm  be  avoided  by  drawing  the  contour  C  slightly  above  the  point 
a~  2 . 

Since  the  behavior  of  the  integrand  of  equation  (61)  has  been  established,  we  apply  the  meth- 
•  of  steepest  descent  (or  saddle  point  method)  to  approximate  the  integral  solution.  This  method 
discussed  in  detail  in  Appendix  F.  Basically,  the  technique  involves  locating  the  point  at  which 
t  •  exponential  of  the  integrand  is  stationary  and  then  adjusting  the  path  through  this  point  so 
t  .  the  function  decreases  rapidly  in  magnitude  on  either  side  of  the  stationary  point. 

Consider  an  integral  of  the  following  form, 


J'^a)  eP^a)da 

Comparing  this  integral  with  equation  (61),  we  conclude  that, 
p  =  kR 

f(a)  =  -i  cos  (a  -  aQ) 

GM  ’  e"T  ViT  (1-8ikr1ro,a)  VW^ 


(64) 

(65) 

(66) 
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For  large  values  of  p,  the  magnitude  of  the  exponential  will  be  very  sensitive  to  changes  in  f(a). 
Thus,  proceeding  as  in  Appendix  F,  the  stationary  point  or  saddle  point  of  our  function  is  located 
by  utilizing  equation  (F-5)  as  follows, 


df{«) 

da 


i  sin  (a  -  a0) 


0 


From  this  last  result,  we  find  that  the  stationary  point  is  given  by, 


a  *  a0 

Hence,  the  major  contribution  of  the  reflected  field  comes  from  those  plane  waves  which  travel 
close  to  path  Rj  in  figure  4.  Also,  as  discussed  in  Appendix  F,  the  path  of  most  rapid  decrease  in 
the  real  part  of  f(a)  is  the  path  of  constant  imaginary  part;  therefore,  the  phase  is  constant  along 
this  path.  For  this  problem  then,  we  will  deform  the  contour  of  integration  so  that  we  are  adding 
those  plane  waves  that  are  close  to  the  specularly  reflected  path  and  that  arrive  with  identical 
phase. 

From  equation  (F-6),  the  path  of  steepest  descent  is  given  by, 
f(a)  =  f(aQ)  -  s2 


This  last  result  can  be  written,  employing  equation  (65),  as, 
cos  (a  -  a0)  =  1  -  is2 


(67) 


The  path  described  by  equation  (67)  is  shown  in  figure  6  as  path  S.  As  long  as  the  saddle  point 
aQ  lies  to  the  right  of  the  branch  point  «bp,  the  path  crosses  the  branch  cut  B  twice  so  that  both 
ends  of  S  lie  on  the  proper  branch.  In  this  case,  the  original  contour  C  can  be  deformed  into  S 
without  any  difficulty. 


Physically,  a@p  represents  the  critical  angle  of  plane  wave  reflection  theory.  For  grazing  angles 
iess  than  agp.  the  radical  in  equation  (54)  becomes  imaginary  and  the  magnitude  of  V(a)  is  unity; 
thus,  the  plane  wave  is  totally  reflected.  If  the  grazing  angle  of  the  plane  wave  is  greater  than 
agp,  the  wave  is  only  partially  reflected.  In  the  following  analysis,  the  asymptotic  solution  to 
equation  (61)  will  completely  describe  the  field  only  if  the  angle  shown  in  figure  4  is  greater  than 
the  critical  angle.  As  point  P  in  figure  4  moves  outward,  the  angle  will  pass  through  the  critical 
angle  and,  from  this  point  on,  another  term  will  appear  in  the  solution.  This  term,  a  result  of  the 
branch  line  integral  associated  with  agp,  will  be  considered  later. 


The  asymptotic  solution  to  equation  (61)  can  be  written  by  substituting  equations  (64), 
(65)  and  (66)  into  equations  (F-10),  (F-14)  and  ( F- 1 5).  Neglecting  terms  of  the  order,  1/(kr)2, 
the  reflected  field  is  given  by, 


‘1*2  3 


e~*^R 
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V(a0> 


Q/"(a0)  -  V'  (aQ)  tan  a 


(68) 


where 

kR  »  1. 
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If  we  had  only  considered  the  first  term,  we  would  have  obtained  the  field  resulting  from  a  single 
plane  wave  following  path  2  in  figure  4.  Indeed,  for  large  kR,  the  second  term  of  equation  (68) 
can  be  ignored  and  the  problem  can  be  treated  by  applying  the  plane  wave  reflection  coefficient 
directly  to  the  spherical  wave  incident  at  the  boundary  between  the  media.  The  second  term  of 
the  solution  provides  a  correction  for  the  plane  wave  reflection  coefficient  when  the  incident 
wave  is  spherical.  Finally,  the  validity  of  equation  (68)  depends  not  only  on  kR  being  large,  but 
also  on  V(a)  being  a  slowly  varying  function.  In  the  derivation  of  the  method  of  steepest  descent, 
it  was  assumed  that  the  exponential  factor,  equation  (60),  decreased  in  magnitude  so  rapidly  on 
either  side  of  the  saddle  point  that  the  remainder  of  the  integrand  given  by  equation  (66)  could 
be  adequately  represented  by  the  first  few  terms  of  its  Taylor  series.  As  can  be  shown  from  equa¬ 
tion  (54),  the  reflection  coefficient  V{a)  varies  rapidly  in  the  vicinity  of  the  critical  angle;  hence, 
equation  (68)  is  not  accurate  for  a0  approximately  equal  to  ct0p,  the  critical  angle. 

Lateral  Wave 

When  the  grazing  angle  ot0  of  figure  4  is  less  than  the  critical  angle  agp,  the  path  of  steepest 
descent  no  longer  crosses  the  branch  cut  twice.  Consequently,  a  circuit  around  the  branch  cut 
must  be  added  to  insure  that  the  path  closes  on  the  same  branch  as  the  original  contour.  Figure  7 
illustrates  this  modified  contour.  The  branch  line  integral  introduced  into  the  solution  corresponds 
to  another  form  of  energy  propagation  which  we  will  now  consider. 

Only  the  reflection  coefficient  V(a)  changes  value  in  crossing  the  branch  cut;  thus,  the  ex¬ 
ponential  factor  of  the  branch  line  integrand  is  identical  to  that  in  equation  (61).  By  introducing 
a  modification  to  the  method  of  steepest  descent,  we  can  evaluate  this  branch  line  integral  asymp¬ 
totically.  Brekhovskikh^  outlined  this  modification  and  obtained  a  complete  solution;  however, 
we  will  only  follow  his  work  to  the  point  of  locating  the  path  of  descent.  Once  the  proper  path 
has  been  found,  it  will  be  transformed  back  to  the  7-plane  so  that  the  results  can  be  applied  to 
the  layered  half-space  problem.  Since  the  branch  cut  must  originate  from  the  branch  point,  we  are 
not  able  to  deform  the  path  as  freely  as  in  the  case  of  contour  C,  figure  7.  We  can,  however,  locate 
the  path  leading  away  from  agp  over  which  the  function  decreases  most  rapidly.  Under  these 
circumstances,  most  of  the  contribution  to  the  integral  comes  from  the  immediate  vicinity  of 
<*bp  and  we  can  approximate  the  integrand  by  a  series  similar  to  that  derived  in  Appendix  F.  In 
the  saddle  point  method,  we  locate  the  saddle  of  the  function,  then  descend  along  the  steepest 
path  into  the  valleys  on  either  side  of  the  saddle.  In  the  present  case,  our  starting  point  is  located 
somewhere  on  the  ridge  flank  or  valley  wall  and  we  descend  along  the  steepest  path  into  the  nearest 
valley. 

The  real  part  ot  f(or)  decreases  most  rapidly  along  the  desired  path;  therefore,  as  shown  in 
Appendix  F,  the  imaginary  part  is  constant.  As  a  result,  the  path  of  descent  is  given  by, 

f(a)  =  f(agp)  ~  s  (0<s<°°)  (69) 

We  can  show  that  this  path  leaves  c*BP  with  a  vertical  slope  and  approaches  the  upper  leg  of  S 
(see  figure  7)  asymptotically.  Figure  8  illustrates  the  branch  line  integral  path  B LI  and  the  steepest 
descent  deformation  Ba  after  transformation  to  the  complex  7-plane.  Also  shown  is  the  transformed 
saddle  point  70  with  its  associated  descent  path  S. 


FIGURE  8  -  Path  of  steepest  descent,  S,  and  branch  line  integral,  BLI,  transformed  to  r-plane. 
Also  depicted  is  the  path,  Ba,  for  asymptotic  evaluation  of  the  branch  line  integral 
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Using  aquations  (51)  and  (52),  we  can  express  tha  branch  lina  intagral  ovar  tha  path  BLI 
j:  follows, 


CJ 

,  i  r  4>*i02  "i 

'5  J,  0  h  1  ~ 


Nota  that  in  tha  above  intagral  tha  usual  replacement  of  J0  by  H0^  has  bean  performed.  From 
figures  6  and  7,  wa  nota  that  no  poles  are  crossed  in  tha  process  of  deforming  tha  branch  cut 
so  that  tha  solution  can  be  written  directly  by  changing  tha  path, 


4>g  ■  2bi  J  H0W  (7r)  e'*^  {l  +  d)  f— 

Ba  L 


2  -  b2  fy2 


If  the  asymptotic  form  for  the  Hankel  function  is  used,  the  above  expression  reduces  to. 


-  a*1*  VT  f  lU+dl1  - 
'  ”  {  01 
oa  L_ 


2  -  b2  /322 


y/T  <h 


Since  the  descent  path,  B*  is  initially  vertical  downward  and  only  the  beginning  of  the  path 
contributes  significantly  to  the  integral,  we  can  assume  that  the  real  part  of  y  does  not  change. 
In  other  words,  y  is  given  by, 

7*k2(1-ix)  (71) 

where  x  is  small  in  the  region  of  interest. 

In  view  of  this  last  conclusion,  we  write  the  following  approximations, 

72  ~  l<2^  (1  -  2ix) 

$2  =  |  k 2 2  -  >2  ~  k2  vTix 


=  V  ki2-72  -  ki  *4/1 - -r 


^12  -  b2  —  01 2 
d)  ik2dx 


35 


N  ADC-81 284-30 


After  substituting  these  approximations  into  equation  (70)  and  transforming  the  integral  to  an 
integration  over  xl?,  the  field  can  be  written  as, 


.  ,-ilkjr^,  U  +  d)l  f  ^.-kjrx  d„ 

012V^  { 

The  integration  may  now  be  performed  and  the  final  result  is, 

^>B  .  ~2'^?  e“i[k2r  +  i*  +  d» 


(73) 


where  Py  is  given  by  equations  (72). 

In  order  to  establish  the  physical  significance  of  the  solution,  consider  the  exponential  factor 
of  equation  (73).  By  employing  the  definition  of  the  critical  angle,  equation  (62),  and  the  geom¬ 
etry  of  figure  9,  we  can  write  equation  (73)  in  the  following  form, 


— 2ibk2  ^-ik2  (r  -  R2  cos  a{jp)  ik f  R2 

7y&  8 


(74) 


The  second  exponential  factor  indicates  that  the  energy  associated  with  the  wave  travels  along  the 
path  with  grazing  angle,  agP'  regardless  of  the  range.  Moreover,  as  shown  in  figure  9,  the  two  legs 
of  this  path,  the  incident  and  reflected  segments,  are  separated  by  a  horizontal  portion  equal  to 
the  range  excess  over  the  minimum  range,  R2  cos  Qtgp.  Since  the  wave  number  is  k2  for  this  hori¬ 
zontal  segment,  the  speed  of  propagation  along  the  interface  is  C7  Not  only  is  the  path  of  this 
wave  significantly  different  from  the  direct  and  specularly  reflected  paths  (refer  to  figure  4): 
according  to  equation  (74)  the  amplitude  decreases  as  r*2.  We  observe  that  this  is  a  much  more 
rapid  decrease  than  the  r~1  decrease  associated  with  spherical  spreading.  Furthermore,  the  shape 
of  the  wavefront  is  a  straight  line  in  the  rz  space,  i.e.,  a  conic  section  in  cylindrical  space,  as  is  seen 
by  allowing  the  exponent  of  equation  (73)  to  assume  a  constant  value, 


k2r  +  ki  \l  1 - -  (z  +  d)  =  const. 

f  c22 

The  line  of  constant  phase  (the  wavefront)  is  given  by, 


z  = 


-1 


tan  agp 


+  const. 


In  contrast,  we  note  that,  from  equation  (68),  the  wavefront  of  the  reflected  wave  (and,  incidental¬ 
ly,  the  direct  wave)  is  a  circle, 


R2  =  r2  +  (z  +  d)2  = 


const. 


This,  of  course,  would  be  the  section  of  a  sphere  in  three  dimensions. 
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FIGURE  9.  —  Onset  of  the  lateral  wave  at  range  greater  than  Rmjn  =  R2  cos  «gp. 
Horizontal  displacement,  L.at,  and  wavefront  are  also  shown. 
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Layered  Half-Space  Branch  Line  Integral 

Now  that  we  have  developed  the  concept  of  a  lateral  wave,  we  can  return  to  the  finite  layer 
problem  and  discuss  the  original  branch  line  integral,  equation  (44).  Unfortunately,  the  problem 
is  more  complex  and,  thus,  we  will  be  forced  to  introduce  some  additional  assumptions. 

By  employing  the  asymptotic  form  of  the  Hankel  function,  equation  (44)  can  be  simplified 
as  follows, 


c2 

*B  =  -  V^"e'T  f  ®‘i7rvTlG(01,02)  -G(0l,-02>]  <*T  (75) 

c2 

This  last  result  may  appear  to  be  in  a  form  suitable  for  the  method  of  steepest  descent;  however, 
there  are  additional  exponential  terms  to  be  considered.  These  additional  terms  arise  from  the 
function  G  in  equation  (75).  If  we  replace  the  sine  and  cosine  functions  in  equation  (30),  the  equa¬ 
tion  that  defines  G,  with  complex  exponentials  and  perform  the  indicated  operations  we  can 
write, 


GlfJi.Uji  {  e-01  <*-<».  e-'lMz  +  d) 

_v  [^-ift  (2H  +  z  -  d)  _  e-i/3-i  (2H  +  z  +  d)  _  e~i/3-,  (2H  -  z  -  d)  +  e~i/3n  (2H  -  z  +  d)~| 

(76) 

-V3  [  ]  +  .  . . 

where  V  =  reflection  coefficient  defined  by  equation  (52). 

In  the  same  manner,  we  discover  that  the  expression  for  G (j3-j ,  -  fa)  is  identical  with  equation 
(76)  with  the  exception  that  j3-|  is  replaced  by 

Before  proceeding  further,  it  is  interesting  to  examine  the  nature  of  the  terms  in  equation 
(76).  Each  of  the  terms,  when  combined  with  the  exponential  from  the  asymptotic  form  of  the 
Hankel  function  assumes  the  following  form, 

e-iK  •  6 1 


+V2  |~"e-i0l  (4H  +  z-d)  _  e-i/3-|  (4H  +  z+d)  _  ^  (4H  -  z -d)  +  ^ift  (4H-z  +  d) 


where 

K,  vector  wave  number  =  yr  +  j3-|  k 
ilf,  position  vector  =  rr  ±  (2nH  ±  z  ±  d)  k 

After  integration  over  the  horizontal  wave  number  7,  this  expression  corresponds  to  the  field  of 
a  spherical  wave  (see  Appendix  A)  at  the  point  given  by  the  position  vector  Ai.  Figure  10  illustrates 
the  correspondence  of  each  term  to  a  spherical  wave  at  an  image  receiver  located  so  as  to  produce 
the  required  number  of  surface  and  bottom  reflections.  Thus,  we  have  shown  that  equation  (24), 
the  general  field  solution  in  the  layer,  can  be  thought  of  as  a  superposition  of  spherical  waves 
travelling  along  the  various  direct  and  reflected  paths  from  source  to  observation  point3. 
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Strictly  speaking,  in  order  to  apply  the  method  of  steepest  descent  to  equation  (75),  we 
must  consider  a  separate  path  for  each  term  in  equation  (76)  and  its  counterpart  for  G(0\,  -  02). 
Instead  of  fo'lowing  such  a  difficult  course,  we  introduce  three  assumptions  to  insure  that  the 
problem  is  tractable.  First,  consider  the  leakage  of  waves  into  the  half-space  for  grazing  angles 
greater  than  the  critical  angle  to  be  significant.  This  appears  to  be  a  reasonable  assumption  in  many 
cases  and  allows  the  neglect  of  terms  multiplied  by  powers  of  the  reflection  coefficient  V.  Further, 
we  shall  assume  that  the  range  is  much  greater  than  the  laigest  value  of  2nH  +  i  +  d  in  the  remaining 
terms  so  that  the  variation  of  the  exponential  factor  in  equation  (75)  dominates  the  integrand. 
Under  these  restrictions,  we  shall  consider  the  path  of  steepest  descent  required  by  the  function, 

e~'7r 

In  this  case,  the  descent  path  is  the  Pekeris  branch  line  path  and,  therefore,  the  approximations 
given  by  equations  (72)  can  be  used.  Notice  that  this  path  is  valid  for  either  the  EJP  branch  line 
integral  or  the  Pekeris  integral;  however,  in  order  to  deform  the  EJP  contour  into  the  descent 
path,  all  of  the  poles  in  the  Pekeris  solution  that  are  crossed  must  be  accounted  for  as  residues. 

Finally,  there  is  one  more  restriction  on  the  approximation  we  are  about  to  make.  The  method 
of  steepest  descent,  in  its  present  form,  breaks  down  when  a  pole  is  located  near  the  saddle  point 
or  the  descent  point.  This  circumstance  is  a  possibility  in  the  layered  half-space  problem.  Whenever 
a  mode  approaches  cut-off,  i.e.,  the  transition  between  trapped  and  leaky  waves,  the  pole  associ¬ 
ated  with  the  mode  approaches  the  branch  point.  Cutoff  is  determined  by  the  characteristic  equa¬ 
tion  for  the  poles,  equation  (33),  which  can  be  written, 

tan  0i  H 


At  the  branch  point  7  =  cj/C2,  and,  consequently,  the  right  side  of  the  above  equation  becomes 
infinite.  In  order  for  a  root  to  exist  at  the  branch  point,  the  left  side  must  also  be  infinite  or, 

cuH 

ci 

7  =  co/c2 

This  is  the  condition  for  cutoff  of  the  n1*1  mode.  We  must  not  use  the  approximation  developed 
below  when  equation  (77)  is  satisfied  for  any  positive  integer  value  of  n. 

Now  that  we  have  outlined  the  restrictions  in  our  theory,  we  may  apply  the  approximations 
indicated  in  equations  (72)  to  determine  the  value  of  equation  (75).  In  addition  to  equations  (72), 
we  observe  that, 

01^  cos^  0iH  +  ^2^  sin^  0iH  :^0i2cos2  0iH  (78) 

provided  that  cos0iH=£O,  since  02  is  much  less  than  0i.  This  condition  is  equivalent  to  requiring 
that  no  modes  are  near  cutoff. 
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Substituting  the  above  approximation,  equation  (78),  and  equations  (72)  in  equation  (75), 
we  can  show  that, 


4ibk23sin/31dsinj31z  -ikpr  f  r-  -kpxr  , 

— — — — o - o - ~e  *  /  n/x  e  *  dx 

v/1<2  cos^ /3-|H  *' 


As  in  the  derivation  of  equation  (73),  this  integral  can  be  evaluated;  hence,  the  asymptotic  field 
resulting  from  the  branch  line  integral  is, 

_  2ibk2  sin  ffid  sin  z  e-ik2r  (79) 

r2  j3-|2  cos^  j3-|H 

In  this  last  result,  0-|  is  evaluated  at  y  =  tj/c2- 

As  in  the  case  of  the  lateral  wave  in  the  two  half-space  problem,  this  field  exhibits  a  spreading 
loss  of  l/r^.  Also,  by  examining  the  exponent  of  equation  (79),  we  see  that  the  phase  velocity  of 
this  wave  is  equal  to  the  sound  speed  in  the  half-space  below  the  layer.  Since  we  have  assumed  that 
the  horizontal  range  is  considerably  larger  than  the  water  depth,  the  two  legs  of  the  lateral  wave's 
path  away  from  the  interface  are  so  small  that  their  effects  are  not  evident  in  equation  (79). 

Application  to  a  Shallow  Water  Waveguide 

In  order  to  illustrate  the  behavior  of  the  field  in  a  fluid  layer,  a  specific  example  will  be  con¬ 
sidered.  The  following  values  were  selected  for  the  parameters  in  figure  1  in  an  attempt  to  describe 
a  region  of  shaliow  water  over  a  fluid  sediment  of  moderate  reflectivity: 

b  =  p]/p2  =  0-5 
C-|  =  1500  m/sec 

C2  =  1 600  m/sec 

d  =  20  m 

H  =  50  m 

z  =  40  m 

Below  cutoff  for  the  first  mode  (21.55  Hz)  the  field  was  determined  by  numerical  integration 
of  equation  (24).  For  higher  frequencies,  poles  appear  on  the  real  axis,  representing  trapped  modes, 
and  the  integration  becomes  extremely  difficult.  Hence,  above  cutoff,  the  field  was  evaluated  by 
summing  the  residues  of  the  real  poles  and  numerically  integrating  the  EJP  branch  line  integral 
given  by  equation  (46).  In  this  manner,  it  was  not  necessary  to  locate  the  complex  poles.  Finally, 
the  lateral  wave  field  was  computed  separately  by  numerical  integration  of  the  Pekeris  branch  line 
integral,  equation  (44).  Each  of  the  numerical  integration  computer  programs  was  patterned  after 
the  algorithm  described  by  Bucket3  which  is  a  Guass-quadrature  scheme  with  adaptive  step  size. 

Figures  11,12  and  13  illustrate  the  changes  in  the  field  as  the  source  frequency  is  increased. 
Propagation  loss,  plotted  as  a  function  of  range,  refers  to  the  ratio  of  intensity  (proportional  to 
pressure  squared)  at  the  observation  point  to  the  intensity  one  meter  from  the  source.  In  decibel 
form  this  ratio  is  simply  20  log/«l»/  since  the  intensity  at  one  meter  was  defined  as  unity  by  the 
left  side  of  equation  (A-5). 
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FIGURE  11.  —  Propagation  loss  in  a  shallow  water  waveguide  for  frequencies  up  to  and  slightly 


FIGURE  12.  -  Propagation  loss  for  only  one  trapped  mode 
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FIGURE  13.  -  Propagation  loss  for  frequencies  bounding  cutoff  of  second  mode  (64.66  Hz) 
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In  figure  11,  the  5  Hz  curve  shows  a  rapid  attenuation  of  the  field  with  range.  This  behavior 
is  typical  of  a  waveguide  well  below  cutoff.  As  the  frequency  increases,  the  curves  remain  smooth 
and  indicate  that  the  field  decays  less  rapidly.  Above  cutoff,  that  is,  above  21.55  Hz,  the  appear¬ 
ance  of  the  first  trapped  mode  induces  oscillations  as  shown  by  the  25  Hz  curve.  Referring  to 
figures  12  and  13,  we  observe  that  the  oscillations  retain  a  similar  form  from  25  Hz  to  60  Hz. 
Effects  of  the  second  trapped  mode,  which  appears  above  64.66  Hz,  are  clearly  evident  as  dis¬ 
tortions  in  the  oscillauons  of  the  75  Hz  curve.  Generally,  as  more  trapped  modes  appear,  the 
shape  of  the  curve  becomes  increasingly  complex. 

Another  significant  component  of  the  field  is  the  contribution  of  the  lateral  wave  as  given 
by  the  Pekeris  branch  line  integral.  This  component  is  plotted  in  figures  14  and  15  illustrating  the 
smooth  decay  with  range  at  all  frequencies.  Of  particular  importance  are  the  curves  at  20  Hz  and 
65  Hz  which  are  substantially  higher  in  level  than  any  of  the  other  curves.  Since  the  frequencies 
associated  with  these  curves  correspond  roughly  to  the  cutoff  frequencies  of  the  first  and  second 
modes,  we  conclude  that  the  lateral  wave  is  strongly  excited  when  a  mode  is  near  cutoff.  In  fact, 
at  20  Hz,  the  lateral  wave  is  the  predominant  field  component  (compare  figures  11  and  14)  and 
this  explains  the  relatively  slow  decay  with  range  although,  at  20  Hz,  the  layer  supports  no  trapped 
modes. 

This  behavior  of  the  lateral  wave  when  a  mode  is  near  cutoff  is  reasonable  both  mathematical¬ 
ly  and  physically.  Mathematically,  a  mode  near  cutoff  implies  that  a  pole  is  close  to  the  starting 
point  of  the  branch  line;  therefore,  the  value  of  the  integral  is  large.  Incidentally,  the  proximity 
of  the  pole  to  k2  renders  the  asymptotic  solution,  equation  (79),  invalid.  In  fact,  the  approxi¬ 
mate  integration  considerably  overestimates  the  field  in  this  case.  Physically,  a  pole  near  k2  reflects 
the  existence  of  a  mode  wh-'se  wavefront  normal  is  close  to  the  critical  angle  at  the  interface. 
From  our  analysis  of  the  two  half-space  problem  and  figure  9,  we  observe  that  this  is  the  preferred 
direction  for  the  interface  wave  in  the  upper  medium. 

One  final  note  concerning  the  data  presented  in  figures  1 1  through  15  is  necessary.  Numerical 
integration  of  the  Pekeris  and  EJP  branch  line  integrals  was  performed  by  employing  the  asymp¬ 
totic  form  of  the  Hankel  function.  As  a  consequence,  the  curves  have  been  constructed  so  that 
the  argument,  7r,  of  each  Hankel  function  is  at  least  2m  Below  cutoff,  the  computation  utilized  a 
polynomial  approximation  to  the  Bessel  function  valid  over  the  entire  range  of  the  Bessel  function 
argument. 
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FIGURE  14.  —  Lateral  wave  contribution  showing  excitation  peak  near  cutoff  of  first  mode  (20  Hz  curve) 


FIGURE  15.  -  Lateral  wave  contribution  showing  excitation  peak  near  cutoff  of  second  mode  (65  Hz  curve) 
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CONCLUSIONS 

The  boundary  value  problem  for  a  point  harmonic  source  in  a  layer  of  fluid  over  a  fluid 
half-space  has  been  solved  by  superposition  of  elementary  wave  functions  and  contour  integration 
in  the  complex  wave  number  plane.  Three  different  forms  of  this  solution  have  been  determined. 
Specifically,  these  are:  (1)  the  basic  superposition  integral  given  by  equations  (23)  through  (25);  (2) 
a  finite  sum  of  real  mode  residues,  equation  (35),  and  the  EJP  branch  line  integral  given  by  equation 
(46);  and  (3)  a  finite  sum  of  real  modes,  an  infinite  sum  of  complex  modes  and  tha  Pekeris  branch 
line  integral,  equation  (45).  All  of  these  solutions  are  mathematically  equivalent;  however,  each  has 
its  own  particular  advantages. 

Direct  numerical  integration  of  the  superposition  integral  is  straightforward  provided  that  the 
layer  is  below  cutoff.  In  the  event  that  trapped  modes  exist,  there  will  be  poles  on  the  real  axis  and, 
as  the  frequency  increases,  the  integrand  becomes  increasingly  oscillatu.y.  As  a  result,  the  numerical 
integration  becomes  difficult.  On  the  other  hand,  poles  on  the  real  axis  are  simple  to  locate  by 
means  of  the  characteristic  equation,  equation  (33);  hence,  their  contribution  to  the  field  is  simple 
to  compute.  In  this  case,  the  EJP  integral  must  also  be  evaluated  and  this  evaluation  is  frequently 
difficult  because  the  integrand  can  oscillate  rapidly.  However,  we  observe  that  the  envelops  of  these 
oscillations  tends  to  have  a  finite  number  of  large  but  narrow  peaks  and  these  may  be  individually 
approximated  by  analytical  functions  and  integrated.  Known  as  virtual  modes  or  resonances,  the 
contributions  from  these  peaks  correspond  in  behavior  to  the  lower  order  leaky  modes  18,19,20. 
The  concept  of  virtual  modes  represents  an  approximation  to  the  EJP  integral  and  yields  a  physical 
interpretation  in  terms  of  leaky  modes  of  propagation. 

In  comparison  with  the  EJP  integral,  the  Pekeris  branch  line  integral  is  simple  to  evaluate 
numerically  12.  An  increase  in  complexity  arises,  however,  as  complex  eigenvalues  must  be  located 
in  order  to  complete  the  solution.  Further,  the  fact  that  each  of  these  complex  modes  becomes 
infinite  with  infinite  depth  is  difficult  to  understand  although  the  field  is  computationally  correct 
10.  Leaky  modes  which  behave  improperly  for  one  coordinate  are  adopted  in  other  fields  such  as 
plate  vibration  21  and  electromagnetic  radiation  22  and  the  series  frequently  converges  rapidly.  In 
addition,  the  Pekeris  branch  line  integral  provides  a  distinct  descriptive  device  for  the  lateral  wave  — 
a  specific  physical  mechanism. 

In  the  course  of  this  investigation,  three  distinct  types  of  propagating  waves  have  been  isolated. 
First,  there  are  a  finite,  sometimes  zero,  number  of  trapped  modes.  These  modes  are,  at  long  range, 
traveling  cylindrical  waves  which  transmit  no  energy  into  the  half-space  below  but,  instead,  spread 
cylindrical ly  in  the  layer.  Because  the  leaky  modes  and  the  interface  wave  decay  more  rapidly  with 
range,  the  field  in  the  layer  at  long  range  may  be  adequately  represented  by  the  trapped  modes 
alone. 

At  shorter  ranges,  the  next  wave  type,  the  leaky  mode,  becomes  important.  These  modes  are 
such  that  their  wavefront  normjis  strike  the  interface  with  a  grazing  angle  greater  than  the  critical 
angle.  Consequently,  energy  is  transmitted  to  the  half-space.  Indeed,  the  greater  the  angle,  i.e.,  the 
higher  the  mode  number,  the  more  energy  is  transmitted  downward  through  the  interface.  Except 
at  very  short  range,  a  manageable  number  of  these  modes  is  adequate  to  characterize  the  leaky  field 
component. 

Finally,  there  is  a  wave  which  travels  along  the  interface,  continuously  radiating  energy  into 
both  the  layer  and  the  half-space.  Usually  the  contribution  from  this  wave  is  small;  however,  if  a 
mode  is  near  cutoff,  this  lateral  wave  can  be  excited  strongly.  Unfortunately,  the  asymptotic  form 
for  this  wave  field,  equation  (79),  is  not  valid  when  a  mode  is  near  cutoff  and  a  numerical 
evaluation  of  the  Pekeris  branch  line  integral  must  be  performed. 
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Tolstoy  23  has  shown  that  tha  present  modal  is  adaquata  to  pradict  acoustic  propagation  in 
shallow  ocean  water  for  both  sinusoidal  sources  and  impulsive  sources  (explosives).  While  the  fluid 
layer  over  a  fluid  half-space  is  a  simplistic  model,  the  concepts  developed  in  this  paper  can  be 
readily  extended  to  more  complicated  problems.  Addition  of  fluid  layers  increases  the  number  of 
simultaneous  equations  in  the  determination  of  the  constants  of  the  depth  equations  (19).  Variation 
of  the  sound  speed  with  depth  changes  the  form  of  these  depth  equations  and  requires  the  use  of 
special  functions  or  numerical  integration;  however,  the  solution  still  reduces  to  a  residue  sum  plus  a 
branch  line  integral.  Introduction  of  elastic  layers  leads  to  another  set  of  modes  and  another  branch 
line  integral  describing  the  shear  wave  field.  While  the  solution  becomes  difficult  to  implement,  the 
concepts  remain  unchanged. 

Not  only  is  the  present  technique  applicable  to  the  solution  of  acoustic  propagation  in  fluids 
and  solids:  it  can  be  adopted  in  other  field  problems.  Vibration  of  submerged  plates  and 
electromagnetic  wave  propagation  have  already  been  mentioned.  The  present  problem  is  also 
analogous  to  the  problem  of  nucleus  stability  in  quantum  mechanics.  In  fact,  the  differential 
equation,  equation  (9),  is  identical  in  form  to  the  one-dimensional  Schroedinger's  equation  24.  By 
analogy,  the  trapped  modes  correspond  to  the  discrete  bound  states  of  particles  “orbiting"  the 
nucleus,  and  the  leaky  modes  correspond  to  radioactive  states  in  which  particles  escape  from  the 
nucleus. 
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APPENDIX  A 

The  expansion  of  a  spherical  wave  in  terms  of  cylindrical  waves  can  be  used  to  transform  a 
point  radiator  into  a  boundary  condition.  A  result  of  this  transformation  is  that  the  complications 
associated  with  the  solution  of  a  non-homogeneous  differential  equation  are  avoided. 

From  the  separation  of  variables  solution  to  the  wave  equation,  the  elementary  solution,  equation 
(13),  to  the  wave  equation  in  axial  symmetry  was  found  to  be, 

ip  =  Jn  (7r)e'(wt  ±0z)  (A-1) 


where 

ip  =  velocity  potential 

7  *  horizontal  wave  number  _ 

0  =  vertical  wave  number  =  V  -  y2 

The  superposition  integral  can  be  written  with  the  horizontal  wave  number  as  a  parameter.  Thus, 
after  dropping  the  time  dependence  of  the  potential  function,  we  have. 


(A-2) 


where 

A  =  arbitrary  function  of  y 
=  time-independent  potential 

Recall  that  a  spherical  wave  can  be  expressed  as  the  general  solution  of  the  wave  equation 
through  the  following  relation, 

,  _  f (cot  -  kR) 
v  R 


where 


f  =  arbitrary  function  describing  the  waveform 
R,  radial  distance  = 

Following  Stratton^  the  wave  function  corresponding  to  a  harmonic  source  i  cylindrical  coor¬ 
dinates  may  be  written  as, 


Pi(cot  -  kR) 

*  "  B - R - 


(A-3) 


After  removing  the  time  factor,  the  wave  function  in  equation  (A-3)  along  the  plane  z  =  0  assumes 
the  form, 


<t>  (r,  o)  = 


3-ikr 


(A-4) 
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Since  equation  (A-2)  represents  the  general  time  independent  solution  to  the  wave  equation,  we 
must  determine  the  function  A(y)  that  reduces  the  general  solution  to  the  spherical  wave  solution. 
By  substituting  z  *  0  into  equation  (A-2)  and  setting  the  result  equal  to  equation  (A-4),  we  con¬ 
clude  that  the  equation  defining  the  function  A(7)  is, 

oo 

“ikr  ^ 

—p —  *  J  A(y)  J0(7r)  d7  (A-5) 

o 

In  order  to  solve  this  integral  equation,  we  use  the  Hankel  transform,  which  is  characterized 
by  the  following  transform  pairZ5( 


(A-6) 


(A-7) 


g(7>  =  f  e'ikr  Jo(7r)  dr 
o 

Moreover,  the  Bessel  function  J0  can  be  replaced  by  its  integral  expansion  to  yield^, 

IT  oo 

9W  -  I  I  e-ir(k-7  cos0)  drd0  (A-8) 

-7 r  o 

The  integration  over  r  is  elementary  provided  that  k  possesses  a  negative  imaginary  part.  If  we 
make  this  imaginary  part  arbitrarily  small,  we  obtain. 
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By  introducing  the  transformation  u  ■  e'5  into  equation  (A-9),  the  integral  may  be  expressed  as 
a  contour  integration  in  the  complex  plane. 


9(7) 


J_  /*  du 

*7  J  u 2  -  2uk/7  +  1 


where  C,  the  contour  *  the  unit  circle. 

To  evaluate  the  integral  using  the  residue  theorem,  we  must  first  determine  the  roots  of  the  de¬ 
nominator  which  are  given  by. 


(A-10) 


u2-2uk/7+1  *  (u-uj)(u-u2) 

*  u2  +  u  (u -]  +U2)  +  ujU2 
Comparing  the  constant  terms  of  this  last  result,  we  find, 
u-|U2  =  1 

It  is  apparent  from  this  result  that  one  root  must  lie  inside  the  unit  circle  and  the  other  must  lie 
outside.  Since  only  the  root  inside  the  unit  circle  contributes  to  the  contour  integral,  we  must 
examine  the  magnitude  of  the  roots.  Therefore,  solving  the  quadratic  equation,  (A-10)  for  u, 


u  =  k  ±  >A2  -  72 


Further,  the  residue  at  the  smaller  root  is  equal  to, 

_ _ 1 _  =  -7 

(u  j  -  u2)  2yk2-72 

Finally,  the  integral  in  equation  (A-9)  can  be  written  as, 

9W  *  7'/-^ 

Introducing  this  last  result  into  equation  (A-7),  we  find  that  on  the  plane  z  =  0, 


(A-1 1 ) 


*-  s 


7  d  Q  (7r ) 
iVk2-72 


d7 


Combining  equations  (A-2),  (A-5),  and  this  last  result,  we  find  that  the  velocity  potential  of  a 
spherical  wave  is, 


\h  =  e'Cjjt 


/ 


J0(7r)  e± 


7d7 


(A-12) 


A -4  ' 


,aL'‘ 
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APPENDIX  B 

In  the  process  of  deriving  the  solution  to  the  field  equation,  equation  (24),  in  terms  of  a  sum 
of  residues  and  a  branch  line  integral,  equation  (45),  we  assumed  that  the  integrals  over  contours 
Ci,  C2,  C3  and  C4,  shown  in  figure  2,  vanish  as  the  radius  of  these  arcs  becomes  infinite.  Indeed, 
we  shall  show  that  this  contusion  is  valid  provided  the  source-to-observation-point  range  is  not 
zero;  that  is,  provided  the  observation  point  is  neither  directly  above  nor  directly  below  the  source. 

Employing  equation  (31),  we  can  write  the  required  integrals  in  the  following  form, 

./■  I  H0^  (7^)  7  G(0i,  02>  <*7  (B-1) 

c1234  c1234 

where  C-|234  =  any  of  the  contours,  C-|,  C2,  C3,  C4  shown  in  figure  2. 

For  very  large  values  of  the  horizontal  wave  number  7,  equations  (18)  reduce,  approximately  to, 

01 ,22  ~  -72 


01,2  ~  ±'7 

Observe  that  the  sign  of  the  above  expression  depends  on  the  location  of  the  contour  with  respect 
to  the  branch  cuts  shown  in  figure  2. 

From  the  correspondence  between  the  0  and  7-planes  illustrated  in  figure  3,  we  notice  that: 

a)  along  and  C2, 

01 ,2  =  '7  (B-2) 

b)  along  03, 

01  =  i7 

(B-3) 

02  ■  -hr 

c)  along  C4, 

01,2  =  -'7  (B-4) 

Also,  along  these  circular  arc  contours,  7  is  given  by, 

7  =  Re'0  (B-5) 

where 

R  co 

0  =  angular  position  along  contour. 


B-2 
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r- 

i  i 


As  a  consequence  of  this  last  result,  the  differential  element,  dy,  can  be  written  as, 
d-y  =  iRe'0d0 


(B-6) 


Further,  the  Hankel  function  in  the  integrand  of  equation  (B-1)  can  be  expressed  in  terms  of  the 
asymptotic  approximation,  equation  (29).  Using  equation  (B-5)  for  the  wavenumber,  we  can  then 
write, 


7  H0^)  (yr)  =s  yiT  (RBi0)Vl  e'*  e  -'rR«'* 


(B-7) 


In  order  to  proceed  with  the  analysis,  we  must  consider  each  contour  separately.  Along  the 
contour  Cl,  the  angle  0  ranges  from  -ir  to  -ir/2  and  the  vertical  wave  numbers,  and  fa,  are  given 
by  equation  (B-2).  Substituting  these  values  for  fa  and  fa  into  the  equation  defining  G(0i,  fa), 
equation  (30),  and  replacing  y  with  the  right  side  of  equation  (B-5),  we  can  show  that, 

G(fa.fa)  *  -IjRe^)*1 

Observe  that  terms  that  decrease  exponentially  as  R  increases  to  infinity  have  been  deleted  from 
this  last  equation.  After  this  last  result  is  combined  with  equations  (B-7)  and  (B-6),  the  contour  in 
tegral  given  by  equation  (B-1)  may  be  written, 

n  _ _  -jt/2 

J  =  2  W_fL  J  J  {ei0)Vi  e  R  tr  Sin  0  +  (z-d)  cos  0] 

rs  ' 


e-iR  [r  cos  0  -  (z-d)  sin  0)  j  ^0 


Consequently, 


-tt/2 


/  *  V2S  /  e  R[rsin0  + (z-d)  cos0]  ^0  (B-8) 

Cl  ^  -7T 

Over  the  range  of  0  given  by  the  limits  of  integration,  the  following  inequalities  are  valid, 

-2rR  ( A  +  A  Q  . 
e  \  *  /  >  erR  sm  0 


2R(z-d) 


M 


e  j  ^  eR(z-d)cos0 

Hence,  equation  (B-8)  can  be  reduced  to  the  following  form, 

-ir/2 


/ 

Cl 


V* 


R(z-d-2r) 


/ 


2R0 


(z-d-r) 


e 


d0 


-7T 
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(B-9) 


Observe  that  the  right  side  of  equation  (B-9)  vanishes  as  R  becomes  infinite.  Thus,  the  integral 
along  the  contour  Ci  can  be  neglected  unconditionally. 

Since  G{0\,  02)  is  an  even  function  of  02  (A  consequence  of  equation  (30))  ana  the  sign  of 
0-\  remains  positive  along  C2  and  C3  (see  equations  (B-2)  and  (B-3)),  equation  (B-8)  can  be  ap¬ 
plied  along  these  contours  to  prove  that, 


Introducing  a  change  of  variables,  we  write, 


x  *  R  cos  0 
xQ  ■  R  cos  0O 
dx  =  -R  sin  0  d  0 


(B-10) 


where 


x  =  distance  along  the  real  7  axis 

xQ  *  distance  from  origin  to  either  72  (for  C2)  or  7j  (for  C3) 
0O  *  upper  limit  of  0  for  appropriate  contour. 

Thus,  equation  (B-10)  becomes. 


e  rR  sin  0  e  x(z-d) 


dx 

R  sin  0 


(B-1 1 ) 


In  addition,  the  following  inequalities  are  valid  over  the  range  of  x  indicated  by  the  limits  of  in¬ 
tegration, 


dx  .  dx 
R  sin  0O  R  sin  0 

e  s'n  $o  >  e  rR  sin  0 

Hence,  equation  ( B-1 1 )  can  be  reduced 


/ 

c2,3 


rR  sin  0O 

<  - 

■v/2nrR  sin  0O 


to  the  folloiwng  form, 
*o 

j "  e  x(z’^)  dx 

o 


B-4 
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Integrating  the  above  and  observing  that,  as  R  becomes  infinite  approaches  -ir/2,  we  write, 


<  -£-f£—  f  .1] 

(z-d)  y^2?rrR  1 


(B-12) 


We  can  show  from  this  last  result  that  the  integrals  over  contours  C2  and  C3  vanish  provided  the 
radius  r  is  not  zero. 

According  to  equation  (B-4),  the  sign  of  is  negative  along  C4;  therefore,  the  expression 
for  G(0i,  0 2 ).  >•&.  equation  (30),  has  the  following  form. 


G(0i,  02 )  "  -IfRe'^)-1  e  cos^»e  sin  0 


Substituting  this  last  result  and  equations  (B-6)  and  (B-7)  into  equation  (B-1),  we  can  write  the 
contour  integral  along  C4  as, 


/  »  -ie  '  4  J  (ei <t>)^  e  R  fr  sin  </>  -  (z-d)  cos  <j>] 

C4  *  -tt/2 


.  e  -iR  [r  cos0  +  (z-d)  sin  <t>]  ^ 


Therefore, 


/  <  /  eR[r  sin  0  -  (z-d)  cos  0]  d<J> 

C4  -tt/2 


(B-13) 


Within  the  range  of  </>  given  by  the  above  limits  of  integration,  the  following  inequalities  are  valid, 
e  — '-t  >  erRsin0 


-R(z-a’)  —  (<(>  +i)  D/  ..  , 

e  ff  v  2  ^  e-R(z-d)cos0 

Hence,  equation  (B-13)  reduces  to  the  form. 


/I 


■RlEdl 


Integrating,  we  see  that, 


y  <  ^rR 

C4 


it  e-nUd)  .  e  -rR 
8rR  r  -  z  +  d 


B-5 


N  ADC-81 284-30 


This  last  result  vanishes  unconditionally  as  R  increases  to  infinity.  Thus,  we  have  shown  that  all 
of  the  circular  arc  contours,  Ci,  C2,  C3  and  C4,  in  figure  2  contribute  nothing  to  the  field  provided 
that  the  radius,  r,  is  not  zero. 


✓ 
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APPENDIX  C 

LOCATION  OF  THE  EIGENVALUES 


C-1 
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Selecting  the  sign  of  the  square  root  by  the  same  argument  used  in  deriving  equation  (21),  we  can 
write, 

-'v1,2  =  h,2  (C‘4) 

Substituting  this  last  result  into  equation  (C-1),  we  find  that, 

+  b  tanh  vj  H  =  0  (C-5) 

If  v  is  complex,  equation  (C-5)  can  be  separated  into  real  and  imaginary  parts, 
i^1  r  -  y-j  j  tanh  t>irH  tan  vijH  +  bi»2r  tanh  -  b^i  tan  i>ijH  =  0 
t>1i  +  v\T  tanh  t»irH  tan  jH  +  bt*2r tan  ^*1  i*-*  +  bt»2i  tanh  rH  =  0 


C-2 
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where 

vyr  *  real  part  of  v-\ 
v|j  =  imaginary  part  of  v\ 
i>2r  -  real  part  of  i>2 

V2\  *  imaginary  part  of  ^2 

Eliminating  tan  from  the  above  two  equations,  so  that  we  may  eliminate  the  periodic  change 
in  sign,  we  may  write, 

i^lr  +  b  V2r  tanh  v\r  H 
t»1i  +  b  t>2i  tanh  J>ir  H 

In  each  region  of  the  7-plane  shown  in  figure  3,  the  signs  of  the  components  of  v\  and  vq  possess 
dennite  values.  Using  this  figure  and  equation  (C-1),  the  signs  can  be  determined  and  introduced 
into  equation  (C-6).  As  a  consequence,  we  find  that  complex  eigenvalues  can  not  exist  in  any  region 
of  the  7-plane  except  region  I.  All  the  other  combinations  of  sign  lead  to  a  positive  quantity  on  one 
side  of  equation  (C-6)  and  a  negative  quantity  on  the  other. 

In  view  of  this  analysis,  we  see  that  the  complex  eigenvalues  are  located  in  the  fourth  quad¬ 
rant  of  the  7  plane  to  the  left  of  the  branch  cut  from  0/C2.  Therefore,  the  EJP  cut  possesses  no 
complex  eigenvalues  since  region  I  does  not  exist  for  this  cut. 

To  show  that  an  infinite  number  of  complex  eigenvalues  exist  for  the  Pekeris  cut,  consider 
the  location  of  eigenvalues  for  large  7.  In  this  case,  equation  (C-3)  can  be  written. 


b  V2\  +  *>i  j  tanh  t"| 
b  t»2r  +  v\t tan^  vy 


r_H 

rH_ 


*1,2  -  T 

Substituting  this  last  result  into  equation  (C-5),  we  discover  that, 

1  sinh  27rH  +  i  sin  27jH 

~b  tanh  71 H  cosh  27rH  +  cos  27jH 

where 

7r  =  real  part  of  7 
7j  =  imaginary  part  of  7 

Since  b  is  real, 

sin  27jH  =  0 

This  last  result  is  equivalent  to  the  statement, 


(C-7) 


(C-8) 


where  n  =  large  integer. 


C-3 
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Introducing  this  last  equation  into  equation  (C-7),  we  can  write, 

p 2  _  sinh  2  7rH 
~~  ~  cos  h  2  7rH  ±  1 

By  solving  the  above  equation  for  7r,  we  discover  that  the  real  part  of  these  eigenvalues  asymp¬ 
totically  approaches  the  value, 

*  -  2H  H  t£t  I  <C'9' 

Thus,  equations  (C-9)  and  (C-8)  describe  the  positions  of  these  complex  eigenvalues  for  large  7. 
In  particular,  equation  (C-8)  demonstrates  that  there  are  an  infinite  number  of  these  eigenvalues. 


/ 


C-4 
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APPENDIX  D 

PLANE  WAVE  EXPANSION  OF  A  SPHERICAL  WAVE 
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APPENDIX  D 

Since  the  behavior  of  plane  waves  at  an  interface  between  two  media  is  well  known  we  will 
examine  the  expansion  of  a  spherical  wave  in  terms  of  plane  waves.  The  derivation  given  below 
follows  Srekhovskik.hlS;  however,  in  order  to  facilitate  comparisons  with  material  in  this  paper, 
the  results  are  presented  in  wave  number  space  rather  than  in  terms  of  the  complex  angle  of  in¬ 
cidence. 

The  elementary  harmonic  plane  wave  function  is  given  by, 

yj,  =  e'(^t  -  K  •  <R)  (D-1 ) 

where 

K,  vector  wave  number  =  7xi  +  7yj  +  0k 
<R,  range  vector  =  xi  +  yj  +  zk 
x,  y,  z  =  observation  point  coordinates 
7X»  7y#  0  =  plane  wave  number  coordinates 

In  accordance  with  the  work  of  Appendix  A,  the  spherical  wave  can  be  written  as  a  superposition 
of  these  elementary  functions  along  z  =  0,  thus, 

OO 

“ “  =  ff  A(7x,7y)  e  +  ?yY>  d7xd7y  (D-2) 

-  OO 

This  last  result  is  a  two-dimensional  Fourier  transform;  therefore,  the  inverse  transform  will  yield^ 


7  =  horizontal  wave  number 

ip  =  angular  coordinate  of  7  shown  in  figure  5. 


SSSMji!1  wi.kihi,!  'iiiiMiiijwarrwjiniwiWi 
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Furthermore,  the  element  of  area  in  tne  transformed  space  is  'drdft  Substitution  of  equations  (D-4) 
in  equation  0-3)  changes  the  integral  form  to, 

2ir 

A(7x,7y!  =“2'  J  f  *~'r[ k  -  7  COS  (0  -  (6)]  d0dr 
o  c 

This  last  integral  is  identical  in  form  with  equation  (A-8)  so  that  the  evaluation  procedure  discussed 
in  Appendix  A  may  be  applied.  Thus,  we  fmd, 


A<7v  7y) 


_ -? _ 


Substituting  this  value  of  A  into  equation  0-2),  we  have. 


e-'kr  _-i_ 

r  2  it 


-i(7xx  +  7y  V) 

7^7 


d7x  d7y 


(D-5) 


Writing  this  last  result  in  terms  of  the  wave  number  vector  and  extending  the  investigation  to  in¬ 
clude  the  z  direction,  we  may  represent  the  velocity  potential  as. 


OO 

r  t  e-i;<  • 

d^d7y 

-  OO 


wnere  ft  vertical  wave  number  =  y/k2  -  72  as  illustrated  in  figure  5. 

T,  s  lost  expression  is  the  expansion  of  a  spherical  wave  in  terms  of  plane  waves. 


(D-6) 


Observe  from  the  geometry  of  figure  5  and  the  fact  that  k  is  a  constant,  that,  as  yx  end  7y 
assume  progressively  larger  values,  0  becomes  imaginary.  Waves  having  imaginary  wave  numbers  are 
called  inhomogeneous  waves.  When  we  consider  the  exponential  part  of  these  waves  (from  equation 
(D-1!),  their  nature  becomes  evident. 


ei(wt -7xX -7yY  +  0z)  =  gi(wt-7xX-7yY)  g  _  a2  (2  >  0) 


where 

a  =  positive  real  number 
<5  =  i  a  . 

Hence,  this  wave  propagates  horizontally  and  is  damped  in  an  exponential  manner  in  both  direc¬ 
tions  vertically  away  from  the  plane,  z  =  0. 

Finally,  equation  (D-6)  is  equivalent  to  the  expansion  of  a  spherical  wave  in  terms  of  cylin¬ 
drical  waves  hs  given  by  equation  (A-12).  This  correspondence  can  be  demonstrated  by  applying 
the  transformation  defined  by  equations  (D-4)  to  equation  (D-6).  As  7X  and  7y  are  altered  from 
-=>o  :c  00,  7  assumes  values  from  zero  to  infinity  and  <p  varies  from  zero  to  2tt.  In  addition,  the 


l  if;. 


D-3 
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differential  element  equivalent  to  d7xd7y  is  7d7d0.  If  we  apply  this  transformation  to  equation 
(D-6),  we  can  write, 


*  = 


-ieifaJt 

2  it 


IS 


e-i7r(sin  0  sin  Q  +  cos  0  cos  6)  e±i/3z 
- 


7d7d0 


OO 


O 


2tt 


y  e-i7rcos(0-0)  d0 


d7 


The  quantity  in  the  brackets  of  the  above  integral  is  equivalent  to  an  integral  form  of  the  Bessel 
function**,  J0  <7r),  therefore. 


=  pi  cut 


\p  =  e 


(7r) 

ip 


ydy 


o 

This  last  result  is  identical  with  equation  (A-12). 


(D-7) 
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APPENDIX  E 

BOUNDARY  VALUE  PROBLEM  FOR  TWO  FLUID  HALF-SPACES 
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APPENDIX  E 


Determination  of  the  field  induced  in  a  fluid  by  a  point  source  ubove  a  fluid  half-space  re¬ 
quires  a  procedure  analogous  to  the  technique  employed  in  the  body  of  this  paper.  Figure  4  illus¬ 
trates  the  physical  problem  in  which  the  tower  half-space  is  a  homogeneous  fluid  of  sound  speed, 
c2,  and  density,  p2.  Depth  functions  corresponding  to  equations  (19)  can  be  written. 


Zi  *  C1  e',<31z 

(d  <  Z  <  ew) 

Z2  -  C2e"'^z  +  C3  e$1  2 

(0  <  z  <  d) 

(E-1) 

Z3  -  C4ei<322 

(Z  <  0) 

Boundary  conditions  applicable  to  this  problem  are:  (a)  continuity  of  pressure  at  z  *  0;  (b) 
continuity  of  particle  velocity  at  z  »  0;  (c)  continuity  of  pressure  at  the  source  depth;  and  (d)  a 
discontinuity  of  particle  velocity  at  the  source  depth  representing  the  point  source.  In  mathematical 
terms,  these  conditions  are, 


(a)  P!  Z2  =  pi  Z3 

jb  dZ2  _  dZ3 

'b>  dz  dz 

(c)  Zt  =  Z2 

(d)  — — ^  =  27 

dz  dz  ^ 


at  z  =  0 
at  z  =  0 
at  z  =  d 

at  z  *  d 


From  these  conditions,  we  can  determine  the  constants  in  equations  (E-1)  ard  write  the  depth 
function  in  the  region  below  the  source  and  above  the  interface. 


22  = 


-i/3-j  (z  +  d)  Pi  ~  -ift  (d 

Pi  +  b02 


\0  <  z  <  d) 


where  b,  density  ratio  =  pi/p2 


Finally,  by  employing  a  superposition  similar  to  that  of  equation  (15),  the  solution  for  the  velocity 
potential  is, 


^2  =  e*4*^ 


J0(7r) 


i/3'1  <d  "  z> 


+  e'wt 


J0(7r) 


,-i.fll  (d 


+  z>  Pi  -  bg2 
^1  +  b/32 
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APPENDIX  F 

METHOD  OF  STEEPEST  DESCEN  T 
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APPENDIX  F 

The  method  of  steepest  descent,  or  saddle  point  method,  is  a  technique. for  approximating 
complex  contour  integrals.  If  the  integrand  satisfies  certain  conditions,  a  region  in  the  complex 
plane  can  be  located  that  contributes  most  of  the  integral's  value.  By  considering  the  poles  and 
branches  of  the  function,  the  original  contour  can  be  deformed  so  that  it  passes  through  this  region. 
Furthermore,  the  path  is  selected  so  that  the  function  decreases  as  quickly  as  possible  on  each  side 
of  this  region.  By  requiring  that  a  parameter  in  the  integrand  be  large,  the  descent  in  either  di¬ 
rection  becomes  steep  and  only  the  behavior  of  the  function  near  a  single  point  need  be  considered. 

Hence,  the  problem  is  one  of  locating  the  saddle  point  where  the  peak  value  of  the  integrand 
occurs,  finding  the  path  of  steepest  descent  away  from  the  saddle  point,  deforming  the  original 
contour  into  the  descent  path  and  writing  a  series  approximation  for  the  integral.3.16,26,27  This 
appendix  foiiows  the  method  of  Brekhovskikh^  although  a  more  direct  technique  for  obtaining 
the  asymptotic  series  is  introduced. 

The  method  of  steepest  descent  may  be  applied  to  integrals  of  the  form, 


G  (a)  ePf<«)  da  (F-1) 


where 

c  =  contour  in  the  a-plane 
G,  f  *  analytic  functions 

I  ■  complex  integral 

a  =  independent  complex  variable 

p  =  positive  real  number 

We  consider  the  evaluation  of  this  integral  for  large  values  of  p. 

Generally,  f  is  a  complex  function  which  can  be  written  as  the  sum  of  a  real  and  an  imaginary 

part, 


f(a)  =  f  i  (a)  +  i  f2<a) 

Furthermore, 

epf  =  (F-2) 

Thus,  fi  controls  the  magnitude  of  the  exponential  while  f2  defines  its  phase.  We  must  find  a 
point,  0i0,  at  which  the  function  fi  is  maximum  and  determine  a  contour  for  which  fi  decreases 
most  rapidly  on  either  side  of  o0. 

In  order  to  locate  a0  oi.d  the  path  of  descent,  we  must  consider  the  properties  of  f  (a)  as  an 
analytic  function  of  '•crnplex  variable  cv.  Since  the  function  is  analytic,  the  Cauchy-Riernann 
conditions^  apply, 


F-2 
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3f 1  df2 
9x  *  9y 


A  stationary  point  of  f  i  is  obtained  if, 
3f1  _  9f1  _  « 


Equations  (F-3)  imply  that  the  above  conditions  are  also  valid  for  the  function  f2>  If  we  assume 
that  fj  is  a  maximum  with  respect  to  x,  then, 

9^f*  9^fn 

— - =-<  0 

9x^  9xdy 

Differentiating  the  second  of  equations  (F-3)  with  respect  to  y,  we  find. 


This  last  result  demonstrates  that  if  f  i  is  a  maximum  with  respect  to  x,  then,  it  is  a  minimum  with 
respect  to  y. 

Hence,  we  conclude  that  for  analytic  functions  a  stationary  point  will  not  correspond  to  an 
absolute  maximum  or  minimum  of  the  function.  If  the  function  f  y  (x,  y)  possesses  a  maximum  in 
one  direction,  we  obtain  a  minimum  at  90°  to  this  direction.  For  this  reason,  a  stationary  point 
of  an  analytic  function  is  called  a  saddle  point  after  the  saddle  surface  to  which  it  corresponds. 
From  the  saddle,  >ve  find  peaks  in  two  opposite  directions  and  valleys  in  the  other  two  perpendicu¬ 
lar  directions.  The  desired  path  of  integration  passes  from  the  saddle  point  through  these  valleys, 
one  on  each  side  of  the  saddle.  This  path  concentrates  the  significant  part  of  the  integration  in  the 
shortest  possible  interval. 

We  recall  that  the  direction  of  most  rapid  change  of  a  function  is  given  by  the  gradient  of  the 
function, 

9fi  dfi 

Vfi  -  -ST  '  +  17  j 

Further,  the  slope  of  this  vector  in  the  x,  y-piane  is  given  by, 

dy  _  dfl/dy 
dx  9fi/9x 

According  to  equation  (F-3),  we  may  write  this  last  result  as, 

dy  _  af2^x 


F-3 
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For  the  path  of  constant  f2,  we  have, 

8fo  3f2 

df2  -  17  +  If  dY  *  0 

This  last  relation  is  identical  to  equation  (F-4).  In  other  words,  the  direction  of  most  rapid  change 
of  f  \  corresponds  to  the  direction  of  constant  f2>  Since  f2  determines  the  phase  of  the  exponential 
in  equation  ( F-2),  the  path  of  steepest  descent  coincides  with  the  path  of  constant  phase. 

Thus,  the  saddle  point  is  determined  by  the  condition, 


Furthermore,  the  required  descent  path  must  oe  a  path  over  which  the  imaginary  part  of  f  remains 
constant  and  the  real  part  decreases  on  both  sides  of  the  saddle  point.  By  introducing  an,  as  yet, 
undetermined  parameter,  we  can  describe  the  behavior  of  f  along  the  descent  path  as  follows, 

f(a)  -  f(ao)  -  s2  (-00  <  s  <  oo)  (F-6) 

where  s  *  real-valued  parameter. 

If  we  can  determine  the  path  in  the  a-plane  along  which  equation  ( F-6)  is  true,  we  have  located 
the  path  of  steepest  descent. 

By  employing  equation  (F-6),  we  can  write  equation  (F*1 )  in  the  following  form, 


/>f(«o) 


J G(a)  e-Ps2 


We  shall  alter  the  variable  of  integration  to  s  by  introducing  a  function  £(s)  which  is  defined  by, 
(•(s)ds  *  G(a)da  ( 

Introducing  equation  (F-8)  into  equation  (F-7),  we  have, 


=  epf(“o>  J  £(s)  e-Ps2  ds 


Since  p  is  large,  only  small  values  of  s  contribute  significantly  to  the  integral.  Therefore,  we  can 
expand  £  in  powers  of  s  and  retain  only  the  first  few  terms  of  such  an  expansion.  In  this  manner, 
equation  (F-9)  becomes. 


I  ^  e 


pf(«0) 


l  7  r  \ 

j  $(o)  J  e-PS2  ds  +y  £"(o)  J  s2e-ps2ds  > 

(  ~  oo  -  oo  f 
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Note  that  terms  with  odd  powers  of  s  when  integrated  yield  zero  since  the  integration  is  between 
limits  symmetric  around  zero.  The  resulting  integrals  have  been  tabulated  in  the  literature  so  that 
the  approximation  to  the  original  integral,  equation  (P-1),  is, 

I  «  ^“o)  VT{  {"(<>))  (P-10) 

To  complete  the  analysis,  we  must  find  an  expression  for  £  in  terms  of  known  functions. 
From  equation  (F-6),  we  may  write, 


Combining  this  last  result  with  equation  (P-8),  we  find, 


(F-11) 


Each  of  the  functions  on  the  right  side  of  equation  (F-11)  can  be  expanded  in  a  Taylor  series, 
thus, 

G(x)  -  G(c)  +  xG'(o)  +~G"(o)  +  ... 

thus,  (F-12) 

f'(x)  «  xf"(o)  +  -y-  f"'(o)  +  y-  f""(0)  +  .  .  . 


where  x  *  a-a0. 

From  equation  (F-6),  we  conclude  that, 

s2  =  --£f''(0)  f"'(o)  f'”'(o) 

Introducing  this  last  expression  and  equations  (F-12)  into  equation  (F-11),  we  have. 


Hence,  the  first  term  of  equation  (F-10)  is  given  by. 


(F-13) 


(F-14) 


Note  that  in  this  last  evaluation,  when  x  =  s  =  0,  then,  a  =  «0.  Differentiating  equation  (F-13) 
twice  with  respect  to  s,  we  obtain  the  second  term  of  equation  (F-10),  i.e.. 


$"<o)  =  2i-(o) 


5(f'")2  +  f,,,' 

12(f")3  4(f")2 


(F-15) 


Observe  that  all  the  quantities  in  equation  (F-15)  must  be  evaluated  at  a  =  ocq. 


F-5/F-6 


